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ABSTRACT 

A  program  involving  both  modeling  and  experiment  has  been  underway  at  Princeton  for  several  years  in  an  effort  to 
explore  the  utility  of  plasmas  and  magnetohydrodynamics  (MHD)  for  aerodynamic  applications.  Early  in  the 
program,  “anomalous”  behavior  of  shocks  on  weakly  ionized  plasmas  has  been  explained  in  terms  of  conventional 
gas  dynamics  with  temperature  gradients.  Then,  theoretical  and  computational  models  have  been  developed  for 
plasma  aerodynamics  and  nonequilibrium  MHD.  The  models  include  a  new  theory  of  nonequilibrium  dissociation 
and  vibrational  relaxation  and  kinetics  of  plasmas  generated  by  electron  beams  and  high-voltage  nanosecond  pulses. 
Aerodynamic  steering  using  plasma  energy  addition  has  been  explored  both  computationally  and  experimentally, 
using  the  newly  developed  microwave-driven  supersonic  plasma  wind  tuimel.  Potential  performance  of  hypersonic 
MHD  devices  with  electron  beam  ionization  has  been  outlined.  In  the  last  year  of  the  program,  plasma  and  MHD 
control  of  hypersonic  flows  and  scramjet  inlets  was  studied.  On-ramp  MHD  device  with  ionization  by  electron 
beams  was  shown  to  be  capable  of  maintaining  the  shock-on-lip  condition  at  Mach  numbers  higher  than  the  design 
one,  while  generating  net  power,  although  at  the  expense  of  some  losses  of  total  pressure.  For  mass  capture  increase 
at  Mach  numbers  lower  than  the  design  one,  a  new  concept  of  “virtual  cowl”  was  proposed  and  tentatively  studied. 
In  preparation  for  experimental  studies  of  MHD  control  of  cold-air  high-speed  flows,  ionization  by  repetitive  high- 
voltage  nanosecond  pulses  was  studied  theoretically,  and  a  plasma  sustained  by  such  pulses  was  made  operational. 
This  plasma  is  to  be  used  for  MHD  studies  at  Mach  4  with  a  6.5  Tesla  superconducting  magnet.  Advanced  optical 
diagnostic  techniques,  including  Filtered  Rayleigh  Scattering,  filtered  Thomson  scattering,  coherent  Rayleigh 
scattering,  and  coherent  Rayleigh-Brillouin  scattering,  have  been  developed  and  applied  to  studies  of  weakly  ionized 
plasmas. 


CHAPTER  1 

SHOCK  PROPAGATION  AND  DISPERSION  IN  GLOW  DISCHARGE  PLASMAS 


Summary 

Spark-generated  shock  waves  were  studied  in  glow  discharges  in  argon  and  argon-nitrogen  mixtures.  UV  Filtered 
Rayleigh  Scattering  was  used  to  measure  radial  profiles  of  gas  temperature,  and  laser  schlieren  method  was  used  to 
measure  shock  arrival  times  and  axial  density  gradients.  Time  accurate,  inviscid,  axisymmetric  CFD  computations 
were  run  and  results  compared  with  the  experiments.  Our  simulation  show  that  changes  in  shock  structure  and 
velocity  in  weakly  ionized  gases  are  explained  by  classical  gas  dynamics,  with  critical  role  of  thermal  and  multi¬ 
dimensional  effects  (transverse  gradients,  shock  curvature,  etc.).  A  direct  proof  of  the  thermal  mechanism  was 
obtained  by  pulsing  the  discharge.  With  sub-millisecond  delay  between  starting  the  discharge  and  shock  launch, 
plasma  parameters  reach  their  steady-state  values,  but  the  temperature  is  still  low,  laser  schlieren  signals  are 
virtually  identical  to  those  without  the  discharge,  differing  dramatically  from  the  signals  in  discharges  with  fully 
established  temperature  profiles. 
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I.  Introduction 


Weakly  ionized  gases  (plasmas)  could  potentially  have  an  impact  on  high-speed  aerodynamics.  Using 
plasmas  for  shock  wave  control,  drag  reduction,  vehicle  steering,  sonic  boom  attenuation,  ignition  of  combustion  in 
engines,  and  MHD  power  extraction  and  enthalpy  by-pass  is  being  actively  discussed.  Of  course,  fundamental  issues 
have  to  be  resolved  prior  to  any  practical  applications. 

Propagation  of  shock  waves  in  weakly  ionized  plasmas  was  studied  in  Russia, and,  later,  in  the 
Shock  acceleration,  broadening,  and  “splitting”  was  observed  in  various  d.c.,  RF,  and  pulsed  discharges,  in  inert  as 
well  as  molecular  gases.  Some  researchers'  "’  attributed  experimental  observations  to  various  plasma-specific 
phenomena,  from  electric  double  layers  to  ion-acoustic  waves  to  new  long-range  interactions  or  microscopic 
stracture  formed  in  gas  discharges.  However,  other  Russian  researchers'^’ found  that  their  data  can  be  adequately 
explained  by  conventional  thermal  phenomena.  Gas  heating,  temperature  inhomogeneity,  and  vibrational  relaxation 
were  theoretically  shown  to  result  in  “anomalous”  behavior  of  shocks  in  gas  discharges.'^'® 

In  Refs.  20  and  21,  analysis  of  various  physical  mechanisms  resulted  in  a  conclusion  that  the  “broadening”, 
“splitting”,  and  acceleration  of  shocks  in  gas  discharges  could  hardly  be  explained  by  plasma-specific  effects.  This 
is  due  ultimately  to  the  very  low  ionization  fraction  in  the  discharges,  which  leads  to  small  fraction  of  the  total  gas 
energy  stored  in  ionized  and  excited  states,  to  very  inefficient  energy  and  momentum  transfer  from  electrons  and 
ions  to  the  bulk  neutral  gas,  and  to  very  long  times  of  these  transfer  processes. 

In  1 997,  two  research  groups  simultaneously  and  independently  showed,  based  on  inviscid  CFD  modeling,^^’ 
that  transverse  and  longitudinal  temperature  and  density  gradients  present  in  gas  discharges  due  to  Joule  heating  can 
be  responsible  for  the  observed  “anomalous”  shock  profiles.  The  inviscid  modeling  of  Ref.  22,  extended  in  Ref.  24, 
emphasized  that  transverse  gradients  result  in  multi-dimensional  nature  of  shock  propagation,  and  that  attempts  to 
interpret  the  phenomena  in  one-dimensional  terms  can  be  misleading.  For  example,  shock  curvature  could  be 
misinterpreted  as  a  “broadened”  or  “split”  ID  structure.  The  inviscid  results  were  confirmed  in  Navier-Stokes 
computations.^^  Additionally,  as  demonstrated  in  Ref.  25,  in  the  case  of  very  weak  shocks  (with  Mach  number  very 
close  to  1)  or  very  long  propagation  distances  (many  meters)  wall  shear  could  also  result  in  an  apparent  “splitting” 
of  laser  schlieren  signals  even  in  a  cold  gas. 

In  this  chapter,  we  describe  a  combined  experimental  and  computational  effort  to  clarify  mechanisms  of  shock 
propagation  in  weakly  ionized  plasmas.  Experiments  were  done  in  well-characterized  positive  column  of  glow 
discharge  in  Ar  and  Ar-N2  mixtures.  Shock  velocity  and  structure  was  measured  using  laser  schlieren  method, 
similar  to  experiments  of  Refs.  12-14.  Since  temperature  and  its  gradients  could  play  a  critical  role  in  shock 
propagation,  radial  temperature  profiles  in  our  experiments  were  measured  with  Filtered  Rayleigh  Scattering. 
Axisymmetric  CFD  modeling  gives  a  very  good  agreement  with  experimental  data,  thus  proving  that  it  is  classical, 
but  multi-dimensional,  gas  dynamics  can  explain  for  the  “anomalous”  effects. 

In  order  to  study  shock  propagation  in  weakly  ionized  plasma  with  very  small  or  absent  temperature  gradients, 
pulsed  discharge  experiments  were  carried  out.  Those  experiments  prove  again  the  thermal  nature  of  shock 
propagation  in  gas  discharges. 

II.  Experiments  on  shock  propagation  in  steady-state  glow  discharges 
II.l.  Experimental  setup  and  procedures 

The  experimental  setup  at  Princeton  University^^’ is  basically  similar  to  that  of  Refs.  12-14.  The  plasmas 
studied  were  created  by  longitudinal  continuos  or  pulsed  glow  discharge.  The  discharge  tube,  3.8  cm  i.d.  and  120  cm 
full  length,  was  made  of  quartz.  Shock  waves  were  generated  by  a  spark  discharge  at  the  one  end  of  the  tube.  The 
spark  circuit  parameters  are:  capacitor  0.2  pF,  voltage  15-20  kV.  The  capacitor  discharge  was  controlled  by  a 
triggering  spark  gap.  The  shock  wave  entered  the  discharge  through  a  narrow  (2  mm  wide  and  0.1  mm  thick)  ring 
electrode  made  of  stainless  steel  foil  mounted  adjacent  to  (almost  flush  with)  the  wall.  The  electrode  nearest  to  the 
spark  was  always  the  anode,  in  order  to  prevent  the  incoming  shock  wave  from  passing  by  the  hot  cathode  and 
through  the  cathode  fall  of  the  discharge.  The  cathode  in  all  pulse  and  most  steady-state  experiments  was  identical  in 
size  and  material  to  the  anode.  In  some  experimental  runs,  where  shock  profiles  deep  inside  the  plasma,  far  from  the 
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cathode,  were  studied,  another  cathode  was  used:  a  25  mm  diameter  hollow  cylinder  with  the  axis  perpendicular  to 
the  tube  axis.  This  cathode  was  used  only  when  it  was  far  enough  downstream  of  the  shock  and  could  not  affect 
shock  propagation.  Both  the  cathode  and  anode  could  be  placed  at  different  location  in  the  tube,  so  that  the  distance 
between  the  spark  and  the  anode  could  be  20  or  40  cm,  and  the  length  of  discharge  could  be  20,  40,  or  60  cm.  The 
flexibility  in  changing  the  discharge  length  was  used  for  measuring  the  electric  field  in  the  plasma  (see  below). 

Similarly  to  Refs.  12-14,  we  used  a  laser  schlieren  method^*  for  detecting  the  shock  wave  front  (Fig.l).  A 
He-Ne  laser  beam  crossed  the  discharge  horizontally  along  the  diameter  and  was  then  focused  onto  the  plane  of  a 
knife  by  a  lens  (f=50  cm).  The  beam  1/e  diameter  was  about  0.3  mm.  A  photodiode  with  a  response  time  of  20  ns 
measured  the  light  intensity  after  the  knife  edge.  The  signal  was  recorded  by  a  digital  oscilloscope  with  1  kQ  resistor 
at  the  input,  and  a  time  resolution  «0.1ps.  In  some  measurements  a  miniature  pressure  sensor  was  also  used. 

To  measure  the  instantaneous  shock  wave  velocity  we  used  an  approach  similar  to  that  proposed  in  Ref.  29. 
After  the  laser  beam  crossed  the  discharge  it  was  retro-reflected  from  a  small  90°  prism  and  passed  through  the 
discharge  once  more  with  the  spatial  offset  and  then  focused  onto  the  knife  edge.  This  produced  two  closely 
separated  schlieren  signals  of  opposite  signs.  Shock  wave  velocity  was  found  by  dividing  the  distance  between  the 
two  beam  passages  (12  to  21  mm)  by  the  time  gap  between  the  schlieren  signals.  This  method  provided  an  accuracy 
of  absolute  velocity  near  3  per  cent  (»10  m/s)  and  that  of  relative  velocity  about  0.5  per  cent  (2  m/s). 

Experiments  were  performed  in  an  Ar-N2  mixture  (1%  of  N2).  Addition  of  nitrogen  substantially  improved 
the  radial  uniformity  of  the  discharge.  A  pure  argon  discharge  contracted  at  currents  i>30  -  40  mA  (for  pressures 
p>30  Torr).  With  nitrogen  the  discharge  was  not  contracted  at  all  conditions  studied  (p<100  Torr,  i<100  mA).  To 
prevent  electrophoresis  and  accumulation  of  impurities,  the  gas  mixture  was  pumped  through  the  discharge  with  a 
rate  of  300  seem,  which  corresponded  to  5  -  10  cm/s  linear  speed  in  the  tube.  Control  measurements  with  the  gas 
flow  turned  off  for  a  short  time  gave  shock  velocities  and  schlieren  profiles  identical  to  those  with  the  flow  on. 

As  was  already  mentioned,  the  flexibility  in  changing  the  discharge  length  allowed  us  to  measure  an 
electric  field  E  in  a  positive  column  of  the  discharge  (cathode  voltage  drop  Vc  could  also  be  found  in  these 
measurements).  For  p=30 — 50  Torr,  i=30-100  mA  the  E  value  was  from  8  to  14  V/cm.  Assuming  a  Bessel  radial 
profile  for  gas  temperature,  we  calculated  the  reduced  electric  field  E/N  (N  is  the  gas  number  density)  and  electron 
number  density  n^.  Averaged  over  the  tube  cross  section,  these  values  are  E/N=1.2-1.4  Td,  n«=(1.2-4)10'®  cm'^  for 
the  condition  range  mentioned  above. 

11.2.  Measurements  of  temperature  profiles 

We  have  used  Ultraviolet  Filtered  Rayleigh  Scattering  (UV  FRS)^°'  to  measure  temperature  profiles  of 
the  discharges  studied.  A  summary  of  these  measurements  is  presented  here.  A  narrow  linewidth  laser  is  used  to 
illuminate  the  sample  gas  and  the  scattered  light  is  imaged  through  a  narrow  band  absorption  filter  onto  a  detector, 
as  shown  schematically  in  Fig.2.  The  laser  should  be  tunable  so  that  it  may  be  tuned  within  an  absorption  notch,  and 
narrow  linewidth  so  that  the  elastic  background  (due  to  windows,  particulates  etc.)  is  narrow  compared  to  the 
spectral  width  of  an  absorption  notch.  The  amount  of  scattered  light  measured  by  the  detector  depends  on  the 
spectral  overlap  of  the  scattered  light  and  the  filter  absorption  profile.  The  filter  offers  very  nearly  100% 
transmission  outside  of  the  absorption  notches,  yet  is  extremely  optically  thick  (suppression  modeled  as  10®  and 
higher)  within  the  absorption  notches.  The  light  scattered  by  the  atoms  and/or  molecules  in  the  flow  is  thermally 
and  collisionally  broadened,  and  is  therefore  spectrally  wider  than  the  laser  and  background.  When  the  laser  is  tuned 
to  overlap  an  absorption  notch,  the  filter  serves  to  suppress  the  background  while  transmitting  a  large  fraction  of  the 
light  scattered  from  the  flow.  By  appropriate  selection  of  filter  conditions,  and  laser  frequency,  one  may  obtain  very 
strong  background  suppression  (greater  than  10®)  while  transmitting  close  to  half  of  the  scattered  light. 

Depending  on  the  conditions  of  a  specific  measurement,  slightly  different  approaches  may  be  used.  In  cases 
of  constant  (known)  pressure  we  use  a  method  which  precludes  the  need  for  background  subtraction.  The  approach 
is  to  ratio  the  FRS  signal  at  a  temperature  to  be  measured,  Tm,  to  the  signal  at  a  known  temperature,  T^f.  Both 
measurements  are  taken  at  the  same  frequency  within  an  absorption  notch  so  that  any  background  is  strongly 
suppressed  in  both  cases.  At  constant  pressure  there  is  a  one-to-one  correspondence  of  temperature  and  density  so 
that  the  measured  ratio  uniquely  defines  the  unknown  temperature  Tm-  The  signal  ratio  is  primarily  determined  by 
the  density  ratio,  with  a  smaller  opposing  contribution  from  lineshape  effects. 
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Modeling  results  for  the  present  measurements  are  shown  in  Figs.  3  and  4.  Fig.  3  shows  the  modeled  filter 
transmission  profile  and  FRS  signal  levels  as  functions  of  frequency.  The  filter  has  a  length  of  5  cm,  a  mercury 
vapor  pressure  of  0.0030  torr,  and  a  temperature  of  315K.  In  this  spectral  region,  mercury  has  six  absorption 
notches  due  to  hyperfine  splitting  and  different  isotopic  contributions.  Fig.  3  is  a  closer  view  of  the  highest 
frequency  absorption  notch  -  the  one  used  in  the  current  work.  The  nominal  location  of  the  lines  is  at  253.7  nm. 
The  modeled  FRS  signal  levels  are  for  a  scatterer  of  50  Torr  of  argon  at  various  temperatures.  The  Rayleigh 
scattering  cross-section  is  not  affected  by  the  weak  ionization  and  so  the  curves  describe  both  argon  gases  and  the 
plasmas  used  in  these  experiments.  The  signals  are  normalized  to  300  K  (by  the  ideal  gas  law)  so  that  their  relative 
magnitude  scales  as  would  be  measured  by  a  detector.  In  the  absence  of  the  filter  the  signal  is  linear  with  density  (or 
1/T),  while  within  the  absorption  features  there  is  also  a  weaker  effect  from  the  variation  of  scattering  linewidth  with 
temperature.  From  such  modeling  one  may  find  the  relationship  between  the  experimentally  measured  FRS  signal 
ratio. 

Ratio  =  FRS  signal  (T=Tm)  /  FRS  signal  (T=Tfef), 


and  the  unknown  temperature  T^.  Fig.  4  is  such  a  curve.  In  this  case  the  unknown  temperature,  Tm,  corresponds  to 
the  “plasma  on”  temperature  while  the  reference  temperature,  Tfef,  corresponds  to  the  “plasma  off’  temperature  -  in 
this  case  T^f  =308  K. 

The  experimental  configuration  for  FRS  temperature  measurements  is  illustrated  in  Fig.  2.  A  high  power, 
tunable,  narrow  linewidth  Titanium  Sapphire  laser^^  was  used  as  the  excitation  source.  The  frequency-tripled  output 
was  delivered  to  the  plasma  tube  through  several  anti-reflection  coated  beam-shaping  optics.  For  these  experiments 
pulse  energies  of  ~  10  mJ  in  the  ultraviolet  (254  nm)  were  used.  The  beam  passed  through  a  quarter  wave  plate  to 
ensure  the  correct  orientation  of  the  linearly  polarized  beam,  and  through  a  half-wave  plate  to  correct  for  the  slight 
elliptical  polarization  introduced  at  the  tube  windows.  Several  irises  were  used  to  reduce  stray  light  carried  with 
beam.  A  long  focal  length  lens  (-90  cm)  was  used  to  weakly  focus  the  beam  to  a  waist  of  -200  microns  at  the 
sample  volume.  The  sample  volume  was  within  the  discharge  tube  and  was  defined  by  the  laser  beam  waist  and  an 
iris  located  between  the  discharge  tube  and  imaging  lens  (f#  =2.4).  The  geometry  was  such  that  the  sample  volume 
had  a  length  of  about  1mm.  The  scattered  light  was  collected  at  a  scattering  angle  of  90  degrees  with  the  imaging 
lens  and  imaged  through  the  mercury  vapor  filter.  The  filter  had  length  5  cm  and  had  two  temperature  controllers 
yielding  a  vapor  pressure  of  0.003  Torr,  and  a  filter  temperature  of  3 1 5  K.  After  the  mercury  filter,  the  light  passed 
through  a  monochromator,  which  acted  as  a  broad  passband  filter,  while  rejecting  plasma  luminosity  and  other 
colors  of  light  from  the  laser.  Finally,  the  scattered  light  signal  was  collected  with  a  Hamamatsu  R960 
photomultiplier  tube.  After  propagation  through  the  discharge  tube  a  quartz  flat  was  used  to  pick  off  a  fraction  of 
the  beam  for  power  normalization  as  well  as  a  frequency  reference.  The  picked  off  beam  illuminated  a  diffuser,  and 
the  elastically  scattered  light  was  detected  through  a  second  mercury  filter  (frequency  reference),  as  well  as 
measured  by  a  photodiode  (power  reference).  The  latter  two  measurements  used  thin  film  interference  filters 
(CORION  G25-254-F)  for  spectral  filtering.  The  three  signals  were  collected  with  a  Stanford  Research  System 
boxcar  onto  a  computer. 

In  plasma  temperature  profile  measurements,  a  thermocouple  was  used  to  the  measure  the  wall 
temperature.  Two  1-D  translation  stages  were  used  to  translate  the  discharge  tube  horizontally  relative  to  the  laser 
beam.  Radial  profiles  along  the  horizontal  axis  were  found  by  translating  the  discharge  tube.  Profiles  were  obtained 
in  the  argon  plasma  as  well  as  the  argon  plus  1%  nitrogen  mixture,  both  at  a  pressure  of  50  Torr,  at  several  different 
values  of  electric  current. 

For  the  range  of  experimental  conditions  studied,  p=30  -  50  Torr,  i=30  -  100  mA,  the  steady-state 
centerline  temperature  range  was  found  to  be  from  440±30  K  to  830±70  K.  Fig.  5  shows  the  experimental  profile  of 
an  argon  plasma  at  pressure  50  Torr  and  current  20  mA.  Figures  6  and  7  are  experimental  profiles  of  an  argon  +1% 
nitrogen  mixture  at  a  pressure  of  50  Torr,  and  currents  of  10  and  40  mA  respeetively.  The  experimental  points  are  in 
good  agreement  with  computation.^'  We  have  plotted  them  here  with  fitted  Gaussian  curves: 


T{r) 

----  =  —  exp 
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where  Tc  and  T„  are  the  centerline  and  wall  temperature,  respectively;  r  is  the  radial  variable;  and  R  is  the  tube 
radius. 
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III.  Results  of  shock  propagation  studies  in  steady-state  discharges  and  comparison  with  CFD  axisymmetric 
modeling 

In  our  earlier  work“'  we  performed  a  2D  inviscid  modeling  of  shock  propagation  through  a  discharge- 
heated  gas  and  concluded  that  transverse  gradient  and  multi-dimensionality  play  a  cracial  role  in  shock  propagation 
in  glow  discharges,  and  that  multi-peak  laser  schlieren  signals  are  simply  due  to  the  shock  curvature.  For  better 
comparison  with  experiments,  in  this  work  we  developed  an  axisymmetric  version  of  our  code.  The  code  solves  the 
axisymmetric  conservative  form  of  the  Euler  equations  on  unstractured,  triangular,  adaptive  meshes.  The  flux 
balance  and  shock  capturing  were  carried  out  by  extending  the  BGK  solver  developed  and  validated  in  Ref.  33  to 
axisymmetric  problems.  The  grid  adaption  technique  was  identical  to  the  one  used  for  the  two-dimensional  code. 
The  “driven”  gas  was  initially  at  rest,  and  a  shock  wave  of  prescribed  strength  separated  initially  the  “driven”  gas 
from  the  “driver”.  The  state  of  the  “driver”  gas  was  specified  according  to  the  Rankine-Hugoniot  conditions. 

The  subsequent  development  of  the  flow  field  was  computed  by  solving  time-dependent  Euler  equations  in 
conservation  form.  Thus,  no  additional  care  was  required  to  satisfy  the  Rankine-Hugoniot  conditions,  and  the 
correct  shock  propagation  and  jump  conditions  across  shocks  and  other  discontinuities  were  captured  by  the 
solution. 

As  found  in  the  earlier  work,^''  to  model  shock  structure  on  a  10-microsecond  time  scale,  simple 
computations  of  a  “piston-driven”,  or  “step”  shock  would  suffice.  This  is  due  to  the  fact  that  the  initial  spark¬ 
generated  shock  broadens  so  much  in  its  20-40  cm  path  to  the  discharge  entrance  or  a  measurement  point  that  on  the 
5-mm  distance  (or  10  ps  time  interval)  from  the  shock  front  the  shock  is  very  close  to  a  “step”  shock.^'*  This 
simplification,  while  being  quite  accurate,  saved  us  from  running  extremely  CPU-intensive  full  computations  of  the 
spark-generated  pulse  shocks.  The  plasma  region  between  the  infinitely  thin  electrodes  was  considered  uniform 
along  the  tube  axis  (x),  and  to  have  a  symmetric  radial  temperature  profile  described  by  the  analytical  fit  (1)  to  the 
experimentally  measured  profile.  Both  wall  temperature  and  the  gas  temperature  outside  the  discharge  were  set 
equal  to  room  temperature.  Thus,  there  was  a  sharp  thermal  diaphragm  between  the  cold  and  hot  regions.  In  absence 
of  diffusion  (both  physical  and  numerical)  this  diaphragm  will  remain  fixed  with  respect  to  a  quiescent  gas  until  the 
incident  shock  impinges  upon  it.  The  initial  position  of  the  incident  shock  in  the  calculations  reported  here  was  one 
tube  radius  to  the  left  from  the  boundary  between  the  cold  and  the  discharge-heated  gases.  Mach  numbers  of  the 
incident  shock  were  selected  so  as  to  give  the  shock  velocity  in  the  uniform  room-temperature  gas  at  a  given 
location  close  to  that  measured  experimentally  with  the  discharge  off. 

For  comparison  with  laser  schlieren  measurements,  the  density  gradient  integrated  in  x  direction  across  the 
narrow  “laser  beam”  and  averaged  in  the  radial  direction  across  the  tube  was  computed.  The  “laser  beam”  width  in 
the  computations  was  adjusted  once  to  match  the  experimentally  determined  width  of  the  schlieren  signal  with  the 
discharge  off.  The  result  of  this  calibration  is  shown  in  Fig.8.  As  seen  in  the  Figure,  computations  reproduce  the 
shape  of  the  experimental  peak  very  well. 

Note  that  in  the  experiments  time  was  measured  from  the  moment  of  spark  firing,  while  the  computations 
started  with  the  shock  located  at  a  distance  of  one  tube  radius  to  the  left  of  the  discharge  region.  Thus,  computational 
time  is  shifted  with  respect  to  the  experimental  one  by  a  constant.  This  is  reflected  in  Figs.  8—11.  For  example,  in 
Fig.  9,  the  experimental  signals  are  plotted  versus  (?  -700  //s) ,  where  t  is  the  time  from  the  spark  firing,  and  the 

con^utational  results  are  plotted  against  (/-420  //s),  with  the  computational  time  t  counted  from  the  moment 

when  the  shock  was  located  at  a  distance  of  one  tube  radius  to  the  left  of  the  discharge  region.  The  constant  shift  of  t 
obviously  does  not  affect  differential  measurements  or  computations  of  shock  velocity. 

Figures  9-11  show  measured  and  computed  schlieren  signals  in  pure  argon  (Fig.9)  and  Ar+1%N2  (Figs. 
10  and  11).  In  each  figure,  upper  plots  correspond  to  the  discharge-off  case,  while  lower  plots  correspond  to  the 
discharge  current  of  20  mA  (Fig.  9),  and  10  and  40  mA  in  Figs.  10  and  11.  Figs.  10  and  11  differ  in  that  the 
schlieren  measurement  point  in  Fig.  11  is  about  11  cm  downstream  of  that  in  Fig.  10.  Measured  and  computed 
velocities  of  shock  waves  are  shown  near  all  the  signals  in  Figs.  9  -  1 1.  As  seen  in  the  figures,  ,  agreement  between 
the  computations  and  experiments,  in  both  shock  velocities  and  the  two-peak  signal  shapes,  is  excellent.  As 
discussed  in  Ref.  24,  the  two-peak  structure  of  the  laser  schlieren  signals  is  due  to  the  curvature  of  the  shock  front  in 
a  region  with  transverse  temperature  gradient.  This  is  illustrated  in  Fig.  12  (a-d),  where  lines  of  constant  density  are 
shown  in  4  moments  of  time,  as  the  shock  propagated  from  the  room-temperature  gas  into  the  plasma.  Note  that 
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flow  perturbations  upstream  of  the  moving  shock,  such  as  the  distortion  of  the  boundary  between  hot  and  cold 
regions  seen  in  Figs.  12  (c)  and  (d),  do  not  affect  shock  propagation.  The  first  peak  in  the  schlieren  signal  comes 
from  the  portion  of  the  shock  that  propagates  through  the  hot  centerline  region.  The  high  temperature  and  low 
density  in  this  region  result  in  both  higher  shock  speed  and  lower  intensity  of  the  peak  compared  with  those 
corresponding  to  colder  near-wall  regions.  The  portion  of  the  shock  moving  through  the  cold  near-wall  region  lags 
behind  and  produces  a  strong  peak  in  the  signal  due  to  the  high  density  near  the  wall.  (Note  that  in  Fig.  12d,  a  well- 
resolved  Mach  stem  is  discernible  near  the  wall).  The  decrease  in  the  schlieren  signal  following  the  first  peak  is  due 
to  the  curvature  of  the  shock.  Indeed,  the  signal  is  proportional  to  the  axial  component  of  the  density  gradient. 
Therefore,  with  the  decrease  of  the  shock  angle  with  the  tube  axis  from  90°  to  lower  values  (Fig.  12d),  the  signal 
weakens  substantially. 

Although  in  Figs.  9  -  1 1  the  agreement  between  computed  and  experimentally  measured  parameters  is 
excellent,  it  is  not  perfect.  Specifically,  computed  shock  velocities  and  distances  between  two  peaks  in  schlieren 
signals  agree  very  well  with  experimental  values.  On  the  other  hand,  computed  height  of  both  peaks  in  the  signals 
relative  to  the  signal  strength  between  the  peaks  is  higher  than  that  in  experiments.  One  of  the  principal  uncertainties 
in  computations  is  due  to  the  use  of  Gaussian  fit  (1)  to  the  experimentally  measured  temperature  profile,  especially 
since  temperature  measurements  have  their  own  errors.  Next,  neither  experimental  nor  computational  (the  latter  - 
due  to  the  mesh  size)  resolution  is  perfect.  Still  another  factor  is  the  use  of  inviscid  model  that  disregards  viscosity 
and  heat  conduction.  According  to  the  data  of  Ref  34,  kinematic  viscosity  of  Ar  at  50  Torr,  300  K,  is  v=2.1  cmVs, 
and  at  50  Torr,  600  K,  it  is  v=7.1  cmVs.  Thermal  diffusivity  values  are,  of  course,  comparable  with  those  of 
viscosity.^'*  The  total  width  of  the  signal  (that  is,  the  distance  between  the  peaks)  in  Figs.  9  -  1 1  is  about  t=10  ps. 
During  this  time,  viscosity  and  heat  conduction  can  smear  velocity  and  density  profiles  on  a  length  scale  of 
S  =  0.1 -0.2  mm  .  Since  the  width  of  each  peak,  converted  into  the  length  scale,  is  about  0.5  mm,  and  the 

distance  between  the  peaks  -  about  5  mm,  viscosity  and  heat  conduction  cannot  change  the  overall  structure  of  the 
schlieren  signals,  which  reaffirms  the  validity  of  inviscid  modeling;  however,  some  smearing  of  the  signals  can 
occur,  helping  to  fill  the  “trough”  between  the  peaks. 

Note  that  since  the  linear  velocity  of  the  gas  downstream  of  the  shock  is  extremely  low,  5-10  cm/s, 
boundary  layer  is  virtually  nonexistent  prior  to  the  shock  arrival.  Behind  the  shock,  boundary  layer  thickness  growth 
with  distance  x  can  be  estimated,  similar  to  the  previous  paragraph,  as  ,  where  t  =  xlu,  and  u  is  the  gas 

velocity  behind  the  shock.  On  the  time  scale  of  1-10  ps  after  the  shock  passage,  the  boundary  layer  will  grow  to  only 
S  =  0.03-0.09  mm  in  our  experimental  conditions.  Thus,  the  boundary  layer  is  very  thin  in  the  vicinity  of  the  shock 
and  does  not  affect  schlieren  signals.  As  to  the  boundary  layer  farther  upstream,  it  has  no  effect  on  the  shock 
propagation. 

For  turther  quantitative  comparison  between  computational  and  experimental  results.  Figures  13  and  14 
show  the  width  of  schlieren  signals  versus  centerline  temperature  (Fig.  13)  and  the  shock  velocity  versus  average 
temperature  in  the  discharge  (Fig.  14),  using  the  same  set  of  shock  propagation  data.  The  average  temperature  was 
determined  in  a  standard  way: 
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Again,  Figure  13  and  14  show  excellent  agreement  between  computations  and  experimental  data.  This  provides  a 
strong  evidence  of  conventional,  thermal  mechanism  of  shock  propagation  in  weakly  ionized  plasmas. 

IV.  Experiments  on  shock  propagation  in  pulsed  discharges 

To  distinguish  between  thermal  and  plasma-specific  mechanisms  of  shock  propagation,  it  would  be 
desirable  to  eliminate  temperature  effects  while  maintaining  plasma  with  the  same  density  of  charged  particles, 
electric  field,  etc.  Unfortunately,  conventional  steady-state  discharge  is  a  system  where  thermal  effects  are  coupled 
with  electric  fields  and  electron  density.  Flowever,  in  pulsed  discharges,  a  relatively  long  time  interval  can  exist 
when  electron,  ion,  and  excited  molecule  densities  are  quite  high  while  the  temperature  is  low.  Indeed,  gas  heating 


6 


and  cooling  occurs  on  time  scale  orders  of  magnitude  longer  than  the  time  scale  for 
ionization/recombination/excitation  processes.^^ 

In  this  work,  the  pulsed  mode  of  the  discharge  was  produced  by  using  a  transistor  switch  in  series  with  the 
discharge.  The  rise  time  for  the  current  pulse  was  20  ps,  and  the  pulse  duration  was  about  0.5  ms.  It  was  found  that 
this  time  was  insufficient  to  get  a  uniform  discharge.  In  fact,  when  the  discharge  was  turned  on,  undesirable 
transitional  processes  (for  example,  discharges  on  the  wall  of  the  tube)  were  observed.  Therefore  a  weak  pilot 
discharge  with  1  mA  current  was  maintained  between  pulses,  which  resulted  in  a  fairly  uniform  volume  pulsed 
discharge.  Fig.  15  shows  the  discharge  pulse  shape  and  the  time  dependence  of  discharge  integral  emission  (with  no 
shock  wave).  Clearly,  at  times  near  to  the  middle  of  the  pulse  the  emission  reaches  its  steady-state  value,  similar  to 
that  of  the  continuous  discharge.  The  initial  peak  of  intensity  is  a  result  of  the  higher  electric  field  arising  in  the 
discharge  immediately  after  the  transistor  switch  is  opened  (see  below). 

As  shown  in  Figs.  16  and  17,  for  the  pilot  discharge  a  small  acceleration,  accompanied  by  some  widening 
and  weakening  of  the  signal  can  be  noticed,  but  the  changes  are  very  small  compared  to  the  higher-current 
continuous  discharge.  This  is  no  surprise  since  the  electron  number  density  in  the  pilot  discharge  is  »10®  cm’^  only 
and  the  measured  axial  gas  temperature  is  less  than  320  K. 

Much  more  important  is  that  schlieren  signal  obtained  from  the  pulsed  discharge,  as  seen  in  Figs.  16  and 
17,  closely  matches  both  no-discharge  and  pilot-discharge  curves,  and  is  very  unlike  the  signal  from  the  continuous 
discharge.  After  the  transistor  switch  is  turned  on,  the  discharge  current  reaches  its  new  steady-state  value  in  »20  ps. 
This  value  is  almost  the  same  as  in  continuous  discharge.  The  electric  field  strength  E  in  the  pulse  could  be  found  by 
subtracting  cathode  drop  V<.  and  voltage  drop  across  the  ballast  resistor  from  the  power  supply  voltage.  The  E  values 
were  found  to  be  somewhat  larger  than  in  the  continuous  discharge,  but  the  E/N  values  are  almost  identical.  Thus, 
the  electron  number  densities  and  mean  electron  energies  (electron  temperatures)  should  be  very  close  to  those  in  the 
continuos  discharge.  This  is  confirmed  by  the  behavior  of  the  discharge  emission  (Fig.  15).  It  reaches  its  steady-state 
value  that  is  very  close  to  that  of  the  continuous  discharge.  An  initial  overshoot  is  related  to  the  over  voltage  applied 
to  the  discharge  in  the  transition  period  (-20  ps)  from  small  to  large  current.  Higher  E/N  values  in  this  interval  result 
in  efficient  excitation  of  atoms  and  molecules,  including  those  with  radiative  lifetime  -1-10  ps  (e.g.,  intensive  bands 
of  the  l"  positive  system  of  the  N2  molecule). 

The  increase  in  gas  temperature,  AT,  during  the  pulse  for  the  pulse  duration  Xp  may  by  calculated  from  the 
simple  balance  equation  iEip  =  CpmiAT,  where  Cp  is  the  specific  heat,  and  mi  is  a  mass  of  gas  per  imit  length.  Here 
we  ignore  heat  losses  because  the  pulse  duration  is  a  factor  of  100  less  than  the  characteristic  time  for  thermal 
conductivity.  For  the  conditions  of  Fig.  16,  this  equation  gives  AT=0.4  K  -  0.6  K.  Thus,  the  pulsed  discharge  has 
electron  component  parameters  (E/N,  n*,  T^)  the  same  as  in  the  continuous  discharge,  but  the  gas  temperature  is  the 
same  as  in  the  pilot  discharge,  that  is,  close  to  room  temperature.  Comparing  the  three  schlieren  curves  of  Fig.  16 
(and,  similarly,  those  of  Fig.  17)  definitely  shows  that  changing  the  electron  density  by  two  orders  of  magnitude  does 
not  affect  shock  wave  propagation,  while  changing  gas  temperature  from  Tg!«320  K  to  Tg«500-600  K  (from  the  pilot 
or  pulsed  discharge  to  the  continuous  discharge)  affects  the  shock  dramatically.  This  result  may  by  considered  as  a 
strong  evidence  of  thermal  mechanism  of  shock  wave  -  plasma  interaction. 

To  further  illustrate  the  role  of  heating  in  shock  dynamics.  Fig.  18  shows  schlieren  signals  taken  with  time 
delay  between  turning  the  discharge  on  and  launching  the  shock  ranging  from  1  ms  to  200  ms.  The  1 -ms-delay 
signal  is  very  close  to  that  in  the  pilot  discharge,  while  the  200-ms-delay  signal  is  close  to  that  in  the  continuous 
discharge.  It  takes  tens  and  hundreds  of  milliseconds  of  the  discharge  existence  for  the  signal  to  get  close  to  its 
steady-state  shape.  The  same  point  is  also  illustrated  in  Fig.  19,  showing  measured  shock  speed  and  schlieren  signal 
width  approaching  their  d.c.  discharge  values  as  the  discharge  is  allowed  to  exist  longer.  Only  gas  heating  in  the 
discharge  takes  that  long;  ionization  and  excitation  processes  are  many  orders  of  magnitude  faster. 

As  the  discharge  is  tinned  off,  one  would  expect  thermal  effects  to  persist  for  a  long  time,  tens  or  even 
hundreds  of  milliseconds,  when  recombination  of  charged  species  and  quenching  of  excited  states  will  be  long  since 
complete.  Indeed,  this  persistence  is  clearly  seen  in  Figs.  20  and  21,  again  supporting  the  thermal  mechanism  of 
shock-plasma  interaction. 

It  should  be  noted  that  vibrational  temperature,  or  populations  of  vibrational  states,  of  nitrogen  molecules, 
could  be  another  parameter  (apart  from  the  gas  temperature)  that  is  quite  different  in  pulsed  and  continuous 
discharges.  The  vibrational  relaxation  characteristic  time  is  1  -  2  orders  of  magnitude  longer  than  the  pulse  duration, 
so  the  vibrational  temperature  in  the  pulsed  discharge  must  be  almost  the  same  as  in  the  pilot  (1  mA)  discharge. 
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To  clarify  the  role  of  vibrational  relaxation,  shock  wave  propagation  through  pure  argon  discharge  was  also 
studied  at  low  currents  where  contraction  did  not  occur.  In  Fig. 22  schlieren  signals  for  the  two  gases  (Ar  and 
Ar+1%N2)  are  compared.  The  main  features  of  schlieren  signals  for  both  cases  are  similar.  For  a  more  detailed 
comparison,  one  has  to  take  into  account  that  addition  of  even  0.1-1%  of  nitrogen  to  argon  dramatically  changes  the 
electron  drift  velocity  for  the  same  E/N.^®  Another  important  factor  found  in  our  experiments  is  that  the  values  of 
E/N  in  Ar+1%N2  are  substantially  lower  than  those  in  pure  Ar.  Since  E/N  in  the  positive  column  is  determined  by 
the  local  ionization  balance, ionization  and/or  recombination  rates  must  be  affected  by  the  addition  of  nitrogen. 
The  exact  mechanism  of  the  decrease  in  E/N  with  nitrogen  addition  is  not  clear  at  this  time,  and,  being  outside  the 
immediate  goal  of  this  work,  can  be  investigated  in  the  future.  One  possible  explanation  is  that  metastables,  such  as 
^2(^2) ,  and  vibrationally  excited  nitrogen  molecules  enhance  the  ionization.  In  any  case,  the  reduction  in  E/N  is 

an  experimental  fact,  and  it  results  in  lower  gas  temperature  in  Ar+1%N2  even  at  the  same  current  as  in  pure  Ar. 
(See,  for  example.  Figs.  5  and  7).  For  a  meaningful  comparison  of  shock  propagation  data  in  Ar  and  Ar+N2 
mixtures,  at  least  one  of  the  parameters  -  electric  current,  electron  number  density,  or  gas  temperature  -  should  be 
the  same  in  both  cases.  Data  presented  in  Fig.  22  for  Ar+1%N2  mixture  are  shown  at  i=20  mA  (the  same  as  for  pure 
argon),  at  40  mA,  where  gas  temperatures  are  close  for  both  discharges,  and  at  60  mA,  where  electron  densities  are 
almost  equal.  It  is  evident  that  for  the  condition  of  almost-equal  gas  temperatures  (40  mA)  similarity  of  the  patterns 
is  the  best.  This  eliminates  unambiguously  the  vibrational  relaxation  effect  and  supplies  one  more  evidence  of  the 
thermal  mechanism  of  shock  dispersion. 

V.  Conclusions 

Extensive  experimental  data  on  shock  propagation  in  steady-state  glow  discharges  in  argon  and  argon- 
nitrogen  mixtures  were  found  to  be  in  excellent  agreement  with  high-accuracy  axisymmetric  CFD  modeling.  This 
provides  compelling  evidence  that  shock  acceleration,  attenuation,  and  “broadening”  in  weakly  ionized  plasmas  can 
be  explained  by  conventional  gas  dynamics,  with  multi-dimensionality  (due  to  transverse  temperature  gradients) 
playing  a  critical  role. 

Experiments  with  pulsed  discharges  allowed  us  to  separate  thermal  effects  from  those  due  to  electric  fields 
and  charged  particles.  It  turned  out  that  gas  heating  is  necessary  for  shock  velocity  and  schlieren  signal  parameters 
to  attain  their  steady-state  values.  This  provides  a  direct  proof  of  thermal  mechanism  of  weakly  ionized  plasma 
effects  on  shock  propagation. 
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Fig.l.  Schematic  of  the  experimental  setup  for  studies  of  shock  propagation  in  glow  discharges 


Fig. 2.  Schematic  diagram  of  experimental  setup  for  UV  FRS  measurements  of  temperature  profiles. 
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Modeled  UV  FRS  signal  levels  for  Argon,  T=300,..,750  K 

Hg  filter;  5cm,  p=0.0030  torr.  T=315  K 


signals  noftnaltzed  to  300  K 
r^otch  used  In  experiments  shown 


Fig.  3.  Modeled  FRS  signal  levels  for  Argon  at  a  range  of  temperatures.  Mercury  has  several  absorption  notches  in 
the  253.7  run  vicinity,  however  Fig.3  shows  the  experimentally  used  (high-frequency)  notch. 


Conversion  Graph:  Signal  Ratio  -  Temperature 
Hg  filter  5  cm,  P=0.0030  torr,  T=315K 


Fig.4.  Look-up  graph  used  to  convert  signal  ratio  found  with  discharge  ON/OFF  to  plasma  temperature. 
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UV  FRS  Temperature  Profile  of  Ar  Plasma,  p=50  torr,  i=20  mA 
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Fig.5.  Radial  temperature  profile  in  Ar  glow  discharge  tube  measured  by  UV  Filtered  Rayleigh  Scattering.  Pressure 
and  electric  current  are  indicated  on  the  plot.  Fitting  curve  -  see  Eqn.  (1). 


UV  FRS  Temperature  Profile  of  Ar+ 1%  N2  Plasma,  p=50  torr,  i=10 

mA 


Fig.6.  Radial  temperature  profile  in  Ar+1%N2  glow  discharge  tube  measured  by  UV  Filtered  Rayleigh  Scattering. 
Pressure  and  electric  current  are  indicated  on  the  plot.  Fitting  curve  -  see  Eqn.  (1). 
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UV  FRS  Temperature  Profile  of  Ar  + 1%  N2  Plasma,  p=50  torr,  i=40 

mA 
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Fig.7.  Radial  temperature  profile  in  Ar+P/oNa  glow  discharge  tube  measured  by  UV  Filtered  Rayleigh  Scattering. 
Pressure  and  electric  current  are  indicated  on  the  plot.  Fitting  curve  -  see  Eqn.  (1). 


t-295  US 

-26  -25  -24  -23  -22 


t-670  ps 


Fig.8.  Calibration  of  the  “laser  beam  width”  for  CFD  modeling  against  laser  schlieren  signal  with  the  discharge  off. 
Note:  experimental  time  t  was  measured  from  the  moment  of  spark  firing,  while  the  computational  time  t  was 
counted  from  the  moment  when  the  shock  was  located  at  a  distance  of  one  tube  radius  to  the  left  of  the  discharge. 
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Fig.9.  Experimentally  measured  and  simulated  laser  schlieren  signals  for  shocks  propagating  in  glow  discharge  in 
pure  Ar  at  50  Torr.  The  discharge  current  is  20  mA.  The  first  of  the  two  laser  beams  is  located  1 8  cm  from  the 
entrance  to  the  discharge,  and  the  spacing  between  the  two  beams  is  3.6  cm.  The  experimentally  measured  and 
computed  shock  velocities  are  indicated  on  the  figure.  Note:  experimental  time  t  was  measured  from  the  moment  of 
spark  firing,  while  the  computational  time  t  was  counted  from  the  moment  when  the  shock  was  located  at  a  distance 
of  one  tube  radius  to  the  left  of  the  discharge. 
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•ig.lO.  Experimentally  measured  and  simulated  laser  schlieren  signals  for  shocks  propagating  in  glow  discharge  in 
y‘+l%N2  at  50  Torr.  The  discharge  current  values  are  10  mA  and  40  mA.  The  first  of  the  two  laser  beams  is  located 
*  cm  from  the  entrance  to  the  discharge,  and  the  spacing  between  the  two  beams  is  2.1  cm.  The  experimentally 
neasured  and  computed  shock  velocities  are  indicated  on  the  figure.  Note:  experimental  time  t  was  measured  from 
he  moment  of  spark  firing,  while  the  computational  time  t  was  counted  from  the  moment  when  the  shock  was 
ocated  at  a  distance  of  one  tube  radius  to  the  left  of  the  discharge. 


Fig.ll.  Experimentally  measured  and  simulated  laser  schlieren  signals  for  shocks  propagating  in  glow  discharge  in 
Ar+1%N2  at  50  Torr.  The  discharge  current  values  are  10  mA  and  40  mA.  The  first  of  the  two  laser  beams  is  located 
20  cm  from  the  entrance  to  the  discharge,  and  the  spacing  between  the  two  beams  is  3.7  cm.  The  experimentally 
measured  and  computed  shock  velocities  are  indicated  on  the  figure.  Note;  experimental  time  t  was  measured  from 
the  moment  of  spark  firing,  while  the  computational  time  t  was  counted  from  the  moment  when  the  shock  was 
located  at  a  distance  of  one  tube  radius  to  the  left  of  the  discharge. 
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Fig.  12.  Computed  constant-density  lines  during  shock  propagation  in  glow  discharge  tube  (pure  Ar,  50  Torr,  20  mA): 
(a)  -  f=21  ps,  prior  to  shock  entry  into  the  discharge;  (b)  —  /=43  ps,  shock  entering  the  discharge;  (c)  —  <=60  ps,  shortly 
after  the  shock  entered  the  discharge;  (d)  -  <=123  ps,  shock  has  penetrated  deep  inside  the  discharge.  IThe 
computational  time  t  was  counted  from  the  moment  when  the  shock  was  located  at  a  distance  of  one  tube  radius  to  the 
left  of  the  discharge 
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Fig.  13.  The  width  of  the  schlieren  signal  in  the  discharge  (A)  minus  the  signal  width  with  the  dicharge  off  (Aq) 
versus  temperature  difference  between  the  axis  and  the  wall,  Ta^is-Tw.  Gas  mixture  Ar+0.16%N2,  pressure  50  Torr;  3 
cm  spark  gap,  15  kV  spark  voltage;  first  laser  beam  is  24.8  cm  inside  the  discharge;  spacing  between  the  two  beams 
is  2.9  cm. 


Fig.  14.  Shock  wave  velocity  versus  average  temperature  [definition  -  see  Eqn.  (2)]  in  the  discharge.  Gas  mixture 
Ar+0. 1  6%N2,  pressure  50  Torr;  3  cm  spark  gap,  1 5  kV  spark  voltage;  first  laser  beam  is  24.8  cm  inside  the 
discharge;  spacing  between  the  two  beams  is  2.9  cm. 
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Fig.l6.  Schlieren  signals  and  shock  velocities  in  glow  discharge  in  Ar+1%N2  at  50  Torr:  the  discharge  off;  38-inA 
continuous  d.c.  discharge,  weak  (1  mA)  pilot  discharge,  and  a  pulsed  discharge. 
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Fig.  17.  Schlieren  signals  and  shock  velocities  in  glow  discharge  in  Ar+1%N2  at  50  Torr:  the  discharge  off;  80-niA 
continuous  d.c.  discharge,  weak  (1  mA)  pilot  discharge,  and  a  pulsed  discharge. 
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Fig.  18.  Evolution  of  the  schlieren  signal  at  a  fixed  location  inside  the  discharge  with  time  delay  between  the 
beginning  of  the  discharge  and  the  shock  launch 
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Fig.l9.  Evolution  of  shock  velocity  and  the  width  of  the  schlieren  signal  af  a  fixed  location  inside  the  discharge  wil 
time  delay  between  the  beginning  of  the  discharge  and  the  shock  launch 
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Fig.22.  Comparison  of  shock  schlieren  signals  at  a  fixed  location  in  Ar  and  Ar+1%N2  discharges  at  50  Torr  and 
various  values  of  electric  current. 
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CHAPTER  2 

THEORETICAL  ANALYSIS  OF  MHD  AND  EHD  EFFECTS  IN  WEAKLY  IONIZED  GAS  FLOWS 


Summary 

This  chapter  focuses  on  theoretical  analysis  of  fundamental  aspects  of  high-speed  flow  control  using 
electric  and  magnetic  fields.  The  principal  challenge  is  that  the  relatively  cold  gas  is  weakly  ionized  in  electric 
discharges  or  by  electron  beams,  with  ionization  fraction  ranging  from  10’*  to  10’*.  The  low  ionization  fraction 
means  that,  although  electrons  and  ions  can  interact  with  electromagnetic  fields,  transfer  of  momentum  and  energy 
to  or  fi-om  the  bulk  neutral  gas  can  be  quite  inefficient.  Analytical  estimates  show  that  even  at  the  highest  values  of 
the  electric  field  that  can  exist  in  cathode  sheaths  of  electric  discharges,  electrohydrodynamic  (EHD),  or  ion  wind, 
effects  in  a  single  discharge  can  be  of  significance  only  in  low-speed  core  flows  or  in  laminar  sublayers  of  high¬ 
speed  flows.  Use  of  multi-element  discharges  would  amplify  the  single-sheath  effect,  so  that  the  cumulative  action 
on  the  flow  can  conceivably  be  made  significant.  However,  Joule  heating  can  overshadow  the  cathode  sheath  ion 
wind  effects.  Theoretical  analysis  of  MHD  flow  control  with  electron  beam  ionization  of  hypersonic  flow  shows  that 
the  MHD  interaction  parameter  is  a  steeply  increasing  function  of  magnetic  field  strength  and  the  flow  velocity. 
However,  constraints  imposed  by  arcing  between  electrode  segments  can  reduce  the  performance  and  make  the 
maximum  interaction  parameter  virtually  independent  of  Mach  number.  Estimates  also  show  that  the  MHD 
interaction  parameter  is  much  higher  near  the  wall  (in  the  boundary  layer)  than  in  the  core  flow,  which  may  have 
implication  for  MHD  boundary  layer  and  transition  control.  The  paper  also  considers  “electrodeless”  MHD  turning 
and  compression  of  high-speed  flows.  Computations  of  a  sample  case  demonstrate  that  the  turning  and  compression 
of  hypersonic  flow  ionized  by  electron  beams  can  be  achieved;  however,  the  effect  is  relatively  modest  due  to  low 
ionization  level. 


Nomenclature 

B  -  magnetic  field 

B.,  -  z-component  and  r-component  of  the  magnetic  field 

B^{z)  -  local  amplitude  of  B-field  on  the  z-axis 
Ttax  ■  maximum  value  of  the  magnetic  field 
e  -  electron  charge 

-  total  energy  (internal  energy,  excluding  the  energy  of  vibrational  mode,  plus  kinetic  energy)  of  the  gas  per  unit 
volume 

-  friction  drag  coefficient 
-  effective  electric  field 

E'  =E  +  uy.B  -  electric  field  in  the  reference  frame  moving  with  the  gas 

£*  -  azimuthal  electric  field  component  in  the  reference  frame  moving  with  the  gas 

^caih  ■  electric  field  strength  at  the  cathode 

E^  -  electric  field  strength  at  the  anode 

E^  -  Hall  field  corresponding  to  the  arcing  threshold 

E^ ,  E^  -  nonequilibrium  and  equilibrium  values  of  vibrational  energy  per  unit  volume 

F_,  -  z-component  and  r-component  of  the  jy.B  force 

f,g-  functions  defined  in  Eq.  (21) 
h  -  flight  altitude 
h  -  Planck’s  constant 
H  -  colutrm  vector  in  Euler  equations 

■  total  enthalpy  entering  the  inlet  per  unit  time 
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j  -  electric  current  density 

-  azimuthal  current  density  component 

_/j  -  current  density  of  the  beam  at  the  injection  point 
k  -  load  factor 

kj  -  rate  coefficient  of  collisional  detachment  of  electrons  from  negative  ions 

k  k  ,  k  ,  k^  ,  k,^,  andi,  -  rate  coefficients  of  collisional  momentum  transfer  between  different  kinds  of 

e+ ”  e—'  en'  h — ”  n— 

particles 

L  -  length  of  MHD  region 

Lg  -  electron  beam  relaxation  length 

M  -  Mach  number 

-  mass  of  an  atom/molecule 

-  positive  ion  mass 
M_  -  negative  ion  mass 
m  -  mass  of  electron 
rk-  mass  flow  rate 

N  -  number  of  consecutive  elements  (multiple  discharges) 

n^,  -  number  densities  of  positive  ions,  negative  ions,  and  electrons 

n  -  number  density  of  gas  molecules 
Nq  -  number  density  of  oxygen  atoms 
p  -  static  gas  pressure 

-  electron  gas  pressure 

p^  -  positive  ion  partial  pressure 
p_  -  negative  ion  partial  pressure 
p^  -  partial  pressure  of  neutral  particles 
Pi-  -  electric  pressure 

<  >  -  total  pressure  averaged  over  the  inlet  cross-section 

Pm.o  ■  free-streara  total  pressure 
q  -  dynamic  pressure 

-  ionization  rate  due  to  electron  beam 

-  total  power  deposited  by  the  e-beam 
q^  -  power  density  deposited  by  the  e-beam 
Qj  -  Joule  dissipation  rate 

-  power  deposited  per  unit  volume  into  vibrational  excitation 


Qyj  -  heating  rate  per  unit  volume  due  to  VT  relaxation 


r  -  coordinate 
-  inlet  radius 


R  -  column  vector  in  Euler  equations 
Rejt/  -  magnetic  Reynolds  number 

Re*  -  Reynolds  number  calculated  at  the  reference  temperature  J* . 


^max  ■  maximal  radius  of  the  numerical  domain 
S  -  MHD  interaction  parameter  (Stuart  number) 

-  MHD  interaction  parameter  with  respect  to  shear  stress  at  the  wall 
T  -  translational-rotational  gas  temperature 
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-  vibrational  temperature 

r*  -  reference  temperature 
U  -  column  vector  in  Euler  equations 
u  -  gas  velocity 

-  freestream  velocity 

M, ,  -  z-component  and  r-component  of  the  gas  velocity 

-  effective  electron-ion  velocity  components  across  magnetic  field  taking  into  account  ion  slip 
u’  -  characteristic  velocity  in  the  near-wall  region  (fnction  velocity) 

V^,  V^,  F  ,  F„  -  convective  velocities  of  electrons,  positive  ions,  negative  ions,  and  neutral  atoms/molecules 
Kr  >  ■  ion  and  electron  drift  velocities 

W-  -  energy  cost  of  ionization,  i.e.  the  loss  of  electron  beam  energy  per  each  newly  generated  electron  in  the  plasma 
X  -  coordinate 

Z  -  column  vector  in  Euler  equations 

z  -  coordinate 

^Met  ■  location 

^max  ■  maximum  length  of  the  numerical  domain 

^EHD  ■  EHD  interaction  parameter 

Z„,  -  total  interaction  parameter 

Zj  -  ratio  of  Joule  and  viscous  dissipation  rates 

Z^  -  EHD  interaction  parameter  for  the  boundary  layer 

a  -  Townsend  ionization  coefficient 

p  -  electron-ion  recombination  rate  coefficient 

-  ion-ion  recombination  rate  coefficient 
Y  -  specific  heat  ratio 

r_-  fluxes  of  electrons,  positive  ions,  and  negative  ions 

-  permittivity  of  free  space 

-  internal  energy  (excluding  the  energy  of  vibrational  mode)  per  unit  mass 

e  -  work  done  on  an  electron  by  the  induced  Faraday  electric  field  during  the  electron’s  lifetime  with  respect  to 
dissociative  recombination  with  ions 

-  initial  energy  of  beam  electrons 

^  -  ratio  of  the  ionization  cost  to  the  work  done  on  an  electron  by  the  induced  Faraday  electric  field  during  the 
electron’s  lifetime 

-  fraction  of  the  Joule  dissipation  rate  spent  on  vibrational  excitation  of  molecules 
6  -  turning  angle  of  the  flow 

-  permeability  of  free  space 
fi  -  dynamic  viscosity 

Vj  -  electron  attachment  frequency  (i.e.  the  number  of  attachments  per  unit  time) 
p  -  gas  density 

-  freestream  density 

a  -  scalar  electrical  conductivity 

Xyj  -  nitrogen  vibrational  relaxation  time 

r  -  shear  stress  at  the  wall 

(p  -  electric  potential 
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-  electron  Hall  parameter 

-  ion  Hall  parameter 


I.  Introduction 


There  is  a  growing  interest  in  using  weakly  ionized  gases  (plasmas)  and  electric  and  magnetic  fields  in 
high-speed  aerodynamics.  Wave  and  viscous  drag  reduction,  thrust  vectoring,  reduction  of  heat  fluxes,  sonic  boom 
mitigation,  boundary  layer  and  turbulent  transition  control,  flow  turning  and  compression,  on-board  power 
generation,  and  scramjet  inlet  control  are  among  plasma  and  MHD  technologies  that  can  potentially  enhance 
performance  and  significantly  change  the  design  of  supersonic  and  hypersonic  vehicles.''^®  Meanwhile,  despite 
many  studies  devoted  to  these  new  technologies,  a  number  of  fundamental  issues  have  not  been  adequately 
addressed.  Any  plasma  created  in  gas  flow  and  interacting  with  electric  and  magnetic  fields  would  result  in  gas 
heating.  This  heating  can  certainly  have  an  effect  on  the  flow  and,  in  some  cases,  can  be  used  advantageously. 
However,  a  more  challenging  issue  is  whether  significant  non-thermal  effects  of  plasma  interaction  with  electric  and 
magnetic  fields  can  be  used  for  high-speed  flow  control. 

In  conventional  magnetohydrodynamics  of  highly  conducting  fluid,  electric  and  magnetic  effects  give  rise 

and  V  —  that  can  be  interpreted  as  gradients  of  electric  and  magnetic  field 


£ 

to  ponderomotive  force  terms  V  — 
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pressures.  These  ponderomotive  forces  are  successfully  utilized  for  plasma  containment  in  fusion  devices  and  also 
play  an  important  role  in  astrophysics.  On  might  hope  that  these  forces  can  also  be  used  for  control  of  high-speed 
flow  of  ionized  air.  However,  the  great  importance  of  ponderomotive  forces  in  fusion  and  astrophysical  plasmas  is 
due  to  the  fact  that  those  plasmas  are  fully,  or  almost  fully,  ionized,  and,  therefore,  are  highly  conductive.  In 
contrast,  high-speed  air  encountered  in  aerodynamics  is  not  naturally  ionized,  even  in  boundary  layers  and  behinds 
shocks  if  the  flight  Mach  number  is  below  about  12,  due  to  the  low  static  temperature.  Therefore,  ionization  has  to 
be  created  artificially,  using  various  electric  discharges  or  high-energy  particle  beams. In  most 
conditions,  the  artificially  created  plasmas  are  weakly  ionized,  with  ionization  fraction  ranging  from  10'*  to  10'*. 
Because  of  the  low  ionization  fraction  and  electrical  conductivity,  interaction  of  the  plasma  with  electromagnetic 
fields  and  transfer  of  momentum  and  energy  to  or  fi'om  the  bulk  neutral  gas  can  be  quite  inefficient. 


In  the  present  paper,  we  derive  analytical  formulas  for  the  so-called  interaction  parameters  that  characterize 
EHD  and  MHD  effects  on  gas  flows.  Numerical  estimates  based  on  these  formulas  are  then  performed  to  determine 
the  trends  and  ranges  of  conditions  where  EHD  and  MHD  effects  can  be  significant. 


II.  Theoretical  analysis  of  MHD  and  EHD  effects  in  weakly  ionized  gas  flows 


2.1.  Basic  equations 

To  analyze  effects  of  electric  and  magnetic  forces  on  weakly  ionized  gas  flows,  we  denote  number  densities 
of  electrons,  positive  and  negative  ions,  and  neutrals  as  n^,  n^,  n_,  and  n ,  velocities  -  as  V^,  V^,  K  ,  and  ,  and 
pressures  -  as  p^,  p^,  p_,  and/;„ .  Rate  coefficients  of  collisional  momentum  transfer  between  different  kinds  of 
particles  will  be  denoted  as  k^_,  k^_,  k^^,  and  k^_ .  Electric  and  magnetic  fields  are  denoted  as  usual:  E  and 

B  .  For  simplicity,  we  will  assume  that  masses  of  both  positive  and  negative  ions  are  equal  to  the  mass  of  a  neutral 
molecule,  and  any  of  those  greatly  exceeds  the  electron  mass  m:  =  M_  □  m  .  Equations  of  motion  for 

four-fluid  mixture  are  then: 
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rf(mF) 

=  -en^  +  F,  X  bJ  - 

'  dt 

- A„«,n^OT(F.  - F, ) - Kn^nm{y^  -  F  )- 

(1) 

d(M„K 

'^-  =  en^{E  +  V^xB)-yp^ 

dt 

(2) 

d(M„t 

^  =  -en  {£-hF 

dt 

- k^n^nm(v_  - F ) - Kn^nM„ (K-V,)- 

-F_«  wM„(f -F„) 

(3) 

=  -^Pr,-k„n^nm(v^-V,) 

dt 

-K,n^nM^{y^-V^)-k„nnM„(v^-V,)  +pAV^  (4) 

where  /i  is  the  dynamic  viscosity.  Adding  these  equations,  we  obtain  a  single-fluid  equation  containing  the  total 
density  p  ,  mass  velocity  u ,  pressure  p  =  p„+  p^+  p^+  p^,  and  electric  current  density  J  =  J : 


p—=-Vp  +  e{n^-n_-nJ)E  +  jx.B+pAu  (5) 

Note  that  the  viscous  terms  in  Eqs.  (4)  and  (5)  are  written  for  the  simplest  case  of  incompressible  flow  with  constant 
viscosity,  because  the  main  focus  is  on  interaction  of  charged  species  with  electric  and  magnetic  fields  and  with 
neutral  molecules  (the  second  and  third  terms  in  the  right-hand  side  of  Eq.  (5)).  However,  the  resulting  Eq.  (5)  can 
be  easily  written  for  the  general  case  of  compressible  flow  with  variable  viscosity:  the  well-known  viscous  terms  of 
Navier-Stokes  equations  should  be  simply  substituted  for  the  last  term  in  the  right-hand  side  of  Eq.  (5). 


2.2,  Electrohydrodynamic  (EHD),  or  “ion  wind”,  flow  control 


Consider  EHD  effects  first.  These  effects,  alternatively  referred  to  as  “ion  wind”,*^  are  represented  by  the 
second  term  in  the  right-hand  side  of  Eq.  (5).  The  strength  of  EHD  effect  is  proportional  to  the  net  space  charge 
density.  Bulk  plasmas  are  known  to  be  quasineutral,  that  is,  in  them,  0  with  very  good  accuracy. 

Substantial  EHD  effects  can  be  expected  only  in  the  presence  of  a  significant  space  charge,  which  is  the  case  in 
cathode  sheaths  or  in  space-charge  regions  of  corona  discharges.  Having  those  space  charge  regions  or  sheaths  in 
mind,  in  what  follows  we  will  neglect  n_  and  in  comparison  with  .  The  change  in  momentum  of  a  gas  element 

moving  along  x  axis  through  a  space  charge  region  is  d(^pu)  =  en^E-dt-en^E-cklu.  Therefore,  the  relative 

strength  of  EHD  effects  can  be  characterized  by  a  dimensionless  EHD  interaction  parameter  equal  to  the  ratio  of 
EHD  push  work  to  the  fluid  momentum  flux: 

_  en^Agr 

^EHD  2 

pu 

where  tk<p  =  ^Edx  is  the  voltage  fall  across  the  space  charge  region. 

An  alternative  way  of  expressing  EHD  effects  is  to  use  the  Poisson  equation,  which  in  a  one-dimensional 
case  can  be  written  as  dEfdx  =  e{n^  —n_-  «,)/£„  .  With  the  Poisson  equation,  the  second  term  in  the  right-hand  side 
of  Eq.  (5)  can  be  transformed  into 

e{n^-n_-n^)E^e^E-j-  =  --  , 

ax  dx\  2  } 
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e  E 

which  can  be  interpreted  as  a  gradient  of  the  “electric  pressure”  ■ .  Therefore,  in  a  one-dimensional  flow 

between  the  anode  and  cathode  of  a  glow  discharge,  the  drop  of  electric  pressure  between  the  anode  and  cathode  is 

where  Ecath  and  are  the  electric  field  strengths  at  the  cathode  and  anode,  respectively.  In  most  cases,  E^  is 
negligible  compared  with  E^^,^  Thus,  the  EHD  interaction  parameter  can  be  expressed  as 


APf  _  ^pEcarh 
pu^  Ipu^ 


(6b) 


Two  principal  examples  of  electric  discharge  systems  that  have  substantial  voltage  drops  across  space  charge 
regions  and  may  be  used  for  EHD  flow  control  are  positive  corona  and  the  cathode  sheath  of  a  glow  discharge. 

In  a  positive  corona  in  atmospheric  air,  can  realistically  reach  about  10®  cm"®,  and  Acp  can  be  from 

several  kilovolts  to  several  tens  of  kilovolts.®®  With  ^(p  =  40  kV ,  the  interaction  parameter  Zehd  will  be  at  least  0.1 

only  for  flows  with  very  low  dynamic  pressure,  3  Pa  or  less.  For  standard  sea-level  air,  this  means  that  the  flow 
velocity  should  be  2  m/s  or  lower. 

In  cathode  sheaths  of  glow  discharges,  ion  density  can  reach  about  lO"  cm"®.®®  Cathode  voltage  fall  in 
normal  (low-current)  glow  discharges  is  about  200-300  V,  and  it  can  reach  about  1  kV  in  abnormal  (high-current) 
glow  discharges.®®  Then,  to  make  >0.1,  the  dynamic  pressure  should  be  less  than  80  Pa,  corresponding  to  the 
maximum  flow  velocity  below  1 1.3  m/s  in  standard  sea-level  air,  or  below  35.8  m/s  in  air  with  1/10  normal  density. 
Note  that  effects  of  surface  glow-like  discharges  with  dielectric  barrier  on  low-speed  flows  have  been 
experimentally  demonstrated.®®’  ®® 

In  principle,  a  stable  discharge  with  electron  density  on  the  order  of  10'®-10'®  cm"®  can  be  maintained 
between  the  cathode  and  the  anode  in  longitudinal  air  flow  preheated  to  2,000-3,000  K  at  a  pressure  of  1  atm.®®  In 
such  systems,  assuming  that  the  temperature  is  2400  K,  the  charge  number  density  is  10’®-10'®  cm"®,  and  the  cathode 
voltage  fall  is  1  kV,  the  EHD  interaction  parameter,  according  to  Eq.  (6a),  can  reach  1  at  dynamic  pressure  of  800- 
8,000  Pa,  corresponding  to  flow  velocities  of  23-72  m/s. 

An  estimate  of  the  maximum  EHD  effect  can  also  be  obtained  with  Eq.  (6b).  If  thermal  emission  of 
electrons  from  the  cathode  is  insignificant,  so  that  electrons  are  emitted  due  to  the  ion  and  photon  bombardment  of 
the  cathode  and  then  multiply  in  the  cathode  sheath  in  Townsend  avalanche  processes,  then  the  cathode  electric  field 
Ec  can  be  very  strong:  » 10“'“'  V  -cm® .®®  Therefore,  for  the  highest  value  of  Zehd  in  Eq.  (6b),  not  only  the 

cathode  sheath  has  to  be  dominated  by  Townsend  processes,  but  the  gas  density  should  be  as  high  as  possible.  If  a 
Townsend  cathode  sheath  could  exist  at  normal  atmospheric  density,  then  the  EHD  interaction  parameter  would 
reach  1  at  dynamic  pressure  of  about  1600  Pa,  corresponding  to  a  flow  velocity  of  about  50  m/s.  At  a  density  of  1/10 
of  standard  sea-level  atmospheric  density,  diffuse  glow  discharges  in  air  with  Townsend  cathode  sheaths  certainly 
can  exist;  in  that  case,  Z  » 1  at  dynamic  pressure  of  about  160  Pa,  corresponding  to  a  flow  velocity  of  about  16  m/s. 

Note  that,  when  assessing  EHD  effects  on  specific  gas  flows,  the  fraction  of  a  gas  actually  passing  through 
the  interaction  region  (i.e.,  through  the  cathode  sheath)  is  as  important  as  the  interaction  parameter  for  that  fraction 
of  the  gas.  Since  at  high  pressures,  where  the  interaction  parameter  can  be  relatively  high,  the  cathode  sheath 
thickness  is  micron-scale  small,  the  EHD  interaction  can  affect  only  a  very  small  portion  of  the  macroscopic  flow. 

Estimates  in  the  previous  paragraphs  referred  to  free  stream  conditions.  However,  from  Eqs.  (6a)  and  (6b)  it 
is  clear  that  maximum  EHD  effects  will  occur  where  the  flow  velocity  is  the  slowest,  that  is,  in  the  boundary  layer. 
For  the  boundary  layer,  the  interaction  parameter  can  be  defined  as  the  ratio  of  electrostatic  and  shear  (friction) 
forces: 


32 


The  wall  shear  stress  can  be  related  to  freestream  conditions  p,  u  with  well-known  correlations:' 


for  laminar  flow 


for  turbulent  flow 


where  Re*  is  the  Reynolds  number  calculated  at  the  reference  temperature  T* .  Therefore, 

z,=^x^  01) 

pu  Cf 

and,  since  « 1  for  high  Reynolds  number  flows,  Z,  greatly  exceed  Z .  For  example,  in  the  range 
Re*  =  10*  -  lO’ ,  the  factor  2/c^  « 1 0’ ,  so  that  ion  wind  effects  in  the  cathode  sheath  of  a  glow  discharge  may  affect 

the  boundary  layer  and  skin  friction  in  standard  sea-level  air  flow  whose  velocity  far  from  the  surface  is  up  to  1 14 
m/s.  In  air  with  0.1  normal  density,  the  velocity  limit  increases  to  360  m/s. 

The  estimates  in  the  preceding  paragraphs,  using  formulas  (6a, b)  and  (7),  assumed  that  a  single  cathode- 
anode  pair  is  used  for  flow  control.  However,  multiple  discharges  can  be  conceivably  used,  so  that  the  cumulative 
effect  of  many  cathode  sheaths  would  affect  the  flow.  A  rough  estimate  of  the  number  of  consecutive  elements,  N, 
required  to  make  the  total  interaction  parameter  Z„,  «  0.1  -1 ,  is  simply  N  =  or  N  =  Z^^JZ^ ,  where  Z^^jj 

and  Zj  are  determined  by  Eqs.  (6a, b)  and  (7). 

Note  that  the  thickness  of  cathode  sheath  of  a  glow  discharge  is  on  the  order  of  a  few  dozen  mean  free 
paths.  Therefore,  only  a  very  thin  portion  of  the  laminar  sublayer  at  the  cathode  would  be  directly  affected.  If  the 
cathode  is  positioned  downstream  of  the  anode,  the  gas  in  this  very  thin  layer  at  the  cathode  would  experience  ion 
wind  forces  directed  both  towards  the  surface  and  in  the  direction  of  the  flow.  What  effect  these  forces  would  have 
on  flow  behavior,  such  as  turbulent  transition  and  heat  transfer,  remains  to  be  studied. 

In  assessing  ion  wind  effects,  it  is  important  to  include  the  inevitable  Joule  heating.  The  ratio  of  Joule  and 
viscous  dissipation  rates  is: 

(.2) 

TU  U 

where  is  the  ion  drift  velocity.  In  the  cathode  sheath,  where  electric  field  is  extremely  strong,  ion  drift  velocity  is 
very  high,  about  (5-10)xl0*  m/s .  Since  u’  is  only  several  meters  a  second,  Eq.  (12)  shows  that  the  heating  effect 

in  the  cathode  sheath  exceeds  the  ion  wind  effect  by  3  orders  of  magnitude.  Thermal  expansion  would  push  the  gas 
from  the  wall,  thus  opposing  the  cathode-directed  ion  wind  force.  The  heating  would  also  increase  viscosity. 
Additionally,  heating  occurs  not  only  in  the  sheath,  but  also  in  the  quasineutral  plasma  region  of  the  discharge,  while 
ion  wind  is  confined  to  the  space  charge  region.  .Thus,  even  if  the  ion  wind  force  can  affect  the  viscous  boimdary 
layer,  its  action  may  be  accompanied  by  the  strong  heating  generated  by  the  same  electric  discharge.  With  regard  to 
the  experiments*’’  **  where  EHD  (ion  wind)  effects  of  surface  glow-like  discharges  with  dielectric  barrier  on  low- 
speed  flows  have  been  demonstrated,  it  would  be  interesting  to  further  investigate  roles  (if  any)  of  viscous  and 
thermal  effects  versus  the  EHD  action  on  the  inviscid  low-speed  core  flow. 


2.3.  Basic  analysis  of  MHD  flow  control 

To  estimate  MHD  effects  in  relatively  cold  hypersonic  flow  and  boundary  layer,  Eq.  (5)  should  be  analyzed 
jointly  with  other  appropriate  gas  dynamic  and  plasma  kinetic  equations.  The  current  density  can  be  related  to  the 
values  of  electric  and  magnetic  field  through  the  generalized  Ohm’s  law  that  is  a  consequence  of  Eqs.  (l)-(4):'‘' 
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(13) 


j  =  aE  +^yxB  +  -^(;xB)x5 
£>  t> 

where  E'  =E  +  uxB  is  the  electric  field  in  the  reference  frame  moving  with  the  gas,  and  and  fi,  are  the 
electron  and  ion  Hall  parameters: 

Q,=-^;Q,=— ^  (14) 

The  second  term  in  the  right-hand  side  of  Eq.  (13)  represents  the  Hall  effect,  and  the  third  term,  nonlinear  with 
respect  to  B,  represents  the  ion  slip.  In  the  analysis,  we  will  disregard  negative  ions  and  also  assume  quasineutrality, 

=  n,  □  n_ . 

The  MHD  interaction  parameter,  also  referred  to  as  the  Stuart  number,  S,  is  the  ratio  of  ponderomotive 
(Ampere)  and  inertia  forces.  With  ion  slip  correction, 

S=  (15) 

(l  +  Qp,)pu 

where  L  is  the  length  of  MHD  region.  When  ion  slip  is  small,  CljEl.  D  1 ,  the  Stuart  number  increases  with  magnetic 
field  as  B^.  However,  at  very  strong  magnetic  fields  and  in  low-density  gases,  when  □  1 ,  the  interaction 
parameter  reaches  its  asymptotic  value  independent  of  B: 

S  -^{Llu)k„^n^  (16) 

The  physical  meaning  of  Eq.  (16)  is  that  the  maximum  ion  momentum  change  occurs  when  they  are  essentially 
stopped  by  the  strong  transverse  magnetic  field,  and  then  momentum  transfer  from  the  ions  to  the  gas  is  limited  by 
the  number  of  collisions  of  a  neutral  molecule  with  ions.  Indeed,  the  right-hand  side  of  Eq.  (16)  is  simply  the  ratio  of 
the  flow  residence  time  in  the  MHD  region  to  the  mean  time  for  a  molecule  to  collide  with  an  ion. 

As  discussed  in  our  earlier  work,'®''*’^^'^‘'’“'^®  at  flight  Mach  numbers  up  to  about  Mach  12,  thermal 
ionization  even  behind  shocks  cannot  provide  an  adequate  electrical  conductivity  for  MHD  flow  control  and  power 
generation  devices.  (Above  Mach  12,  one  might  achieve  adequate  ionization  near  stagnation  zones  and  in  the 
boundary  layers,  but  adequate  ionization  will  not  exist  in  the  inviscid  flow  of  slender  bodies  for  Mach  numbers 
much  higher  than  12.)  The  ionization  must  be  created  and  sustained  artificially,  and  the  need  to  spend  power  on 
ionization  severely  limits  performance  of  hypersonic  MHD  devices. Indeed,  the  work  done  on  an 
electron  by  the  induced  Faraday  electric  field  during  the  electron’s  lifetime  with  respect  to  dissociative 
recombination  with  ions  is: 


mk^Jnn^  mk^Jnn^ 


where  k  is  the  load  factor,  E  -  kuB  is  the  induced  electric  field,  =  - 


is  the  electron  drift  velocity,  and  p 


is  the  dissociative  recombination  rate  constant.  Although  in  the  general  problem  of  MHD  flow  control  one  might 
also  be  interested  in  the  case  where  the  device  is  not  self-powered,  self-powered  MHD  devices  are  obviously 
attractive,  since  they  would  not  require  any  on-board  power  input.  In  this  paper,  we  limit  our  analysis  to  self- 
powered  MHD  devices.  For  efficient  self-powered  operation  of  the  MHD  device,  e  must  be  substantially  larger 
than  the  energy  cost,  ,  of  a  newly  produced  electron.  Specifically,  the  ratio 
mk^JnnJV, 

must  be  limited  to  a  number  less  than  1.  This  obviously  limits  the  electron  density  and  the  conductivity 
He^kuB) 

a  =  =  ^ ^  (19) 

{mk^„n)  pw, 

The  maximum  interaction  parameter  per  unit  length  is  then 
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X - ^  --x^  (20) 

mKnK.P^ 

Hypersonic  MHD  devices  would  operate  downstream  of  one  or  more  oblique  shocks,  and  it  is  convenient 
to  express  the  interaction  parameter  in  terms  of  flight  velocity,  Mach  number,  and  dynamic  pressure.  Assuming  that 
the  MHD  device  is  located  downstream  of  a  single  oblique  shock  with  flow  turning  angle  9 ,  we  express  the  density 
and  velocity  in  terms  of  the  ffeestream  parameters  M,  u^,  with  approximate  formulas  that  can  be  easily  derived 
for  oblique  shocks  with  small  turning  angles: 

P  =  pJf'^  ~  =  -fx — ; 

P  P^  8  _ 

- f - (21) 

'  (y+l)  (y  +  l)Mtan^  +  y(y  +  l)^M^tan^0  +  16 

.  ^  1  AMtznd 

g{M,0)  =  \ - ^  - .  - 

(y  +  l)Mtan0  +  -^(;'  +  l)  M^tan^0  +  16 

With  Eq.  (21),  and  expressing  the  freestream  density  in  terms  of  dynamic  pressure  q  and  velocity,  p^  =  2g/M^  ,  we 
obtain  the  final  expression  for  the  maximum  Stuart  number  per  unit  length: 


L  I,  J  SgjSW^q^  [  J 

MHD  performance  as  expressed  by  Eq.  (22)  is  inversely  proportional  to  the  energy  cost  of  producing  an 
electron,  .  Thus,  minimization  of  W.  is  critical.  As  shown  in  our  earlier  work,*®  '®’^®’^®'^^'^^  DC,  RF,  and 

microwave  methods  of  plasma  generation  have  unacceptably  high  ionization  cost,  If].  0  10“’  eV  .  Electron  beams 
represent  the  most  efficient  nonequilibrium  method  of  ionization,  ’  '  ’  with  fV^  =  34  eV ,  and  the  concept  of  cold- 
air  hypersonic  MHD  devices  with  ionization  by  electron  beams  was  suggested  and  developed  in  our  earlier  work.'®' 
In  what  follows,  we  will  assume  that  ionization  is  done  by  electron  beams,  with  Wg  =  34  eV ,  and  that  ^  is  limited 
to  0.3. 

The  maximum  MHD  interaction  parameter  (22)  increases  very  rapidly  (especially  if  the  ion  slip, 
represented  by  the  last  factor  in  Eq.  (22)  is  not  significant)  with  increasing  magnetic  field  and  flight  speed,  and  with 
decreasing  dynamic  pressure.  The  very  sharp  m’  dependence  is  due  to  both  increase  in  the  Faraday  e.m.f.  with  flow 
velocity  and  to  the  decrease  in  gas  density  at  constant  q.  Figure  1  shows  S/ L  calculated  with  Eq.  (22)  versus  flight 
Mach  number  at  4  different  magnetic  field  strengths  and  2  values  of  flight  dynamic  pressure.  The  load  factor  value 
was  constant  at  ^0.5. 

At  low  gas  density  and  strong  B  field  that  favor  MHD  performance,  Hall  effect  becomes  very  significant.  If 
Hall  current  is  allowed  to  flow,  the  effective  conductivity  is  reduced  by  the  factor  (l  +  Q^J,  and  the  MHD 

performance  decreases.  In  Faraday  MHD  devices,  Hall  current  is  eliminated  by  segmenting  electrodes."'^  However, 
as  the  longitudinal  Hall  electric  field  increases,  arcing  between  the  electrode  segments  can  occur.^^  The  arcing  would 
essentially  result  in  a  continuous-electrode  Faraday  device,  with  dramatic  reduction  in  performance. 

Denoting  the  value  of  Hall  field  corresponding  to  the  arcing  threshold  as  ,  an  additional  constraint  in 
assessing  the  maximum  performance  is: 

£,  =  (1  -  k)n^uB  =  a—  (23) 

ken^P 

This  condition  limits  the  allowed  B  field  at  high  flow  velocity  and  low  gas  density.  With  this  constraint,  the  upper 
limit  for  the  interaction  parameter  per  unit  length  becomes 
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^lin,  _ 

L 


f 


1  + 

V  2(1 


(24) 


This  formula  does  not  explicitly  depend  on  B  field  and  has  a  weak  dependence  on  flight  conditions,  because  at  each 
flight  regime,  B  field  is  adjusted  to  satisfy  Eq.  (23). 

Fig.  1  shows  the  Hall  field-limited  interaction  parameter  of  Eq.  (24)  at  two  values  of  .  While  the 
performance  is  quite  good  at  =  50  kV/m ,  it  is  unacceptably  low  with  E^  =5  kV/m .  Therefore,  determining  the 
threshold  field  for  intersegment  arcing  is  critical  for  assessing  performance  of  hypersonic  MHD  devices. 

The  intersegment  arcing  in  conventional  high-temperature  MHD  generators  with  thermal  ionization  was 
extensively  studied  both  theoretically'*^"'’^  and  experimentally.'’’’'’*  Computational  study  by  Oliver'’^  established  40  V 
as  the  intersegment  potential  drop  corresponding  to  the  onset  of  arcing,  and  experiments'’*  gave  the  range  of  40-100 
V.  With  centimeter-scale  electrode  width  and  intersegment  gap,  this  would  correspond  to  the  longitudinal  field  of  50 
V/cm  or  so.  However,  experimental'”’'’*  and  theoretical'’'’’'’*  studies  of  this  thermally  induced  arcing  found  that  the 
arcing  threshold  increases  with  increasing  magnetic  field  and  decreasing  electrode  surface  temperature. 
Qualitatively,  this  can  be  explained  by  the  reduction  in  electrical  conductivity  in  the  boundary  layer  as  the  wall  gets 

colder,  and  by  the  effective  boundary  layer  conductivity  reduction  as  (l  +  O’ )  ' .  Thus,  even  in  conventional  MHD 


devices  with  thermal  ionization  E^  can  be  substantially  higher  than  50  V/cm. 

More  importantly,  the  operating  conditions  and  the  ionization  mechanism  in  cold-air  hypersonic  MHD 
devices  differ  dramatically  from  those  in  conventional  thermal  systems.  In  on-ramp  MHD  shock  control  and  power 
generation  devices,”’*’  ionization  (n^OlO'*  cm'^,  0.5  mho/m)  is  sustained  by  electron  beams.  The  rate  of 

ionization  by  plasma  electrons  is  very  low,'*"'*  because  the  ratio  of  effective  electric  field  strength  to  the  gas  number 
density,  10'’*  V  ■  cm*  both  in  the  core  flow  and  in  the  boundary  layer,  results  in  low  electron  temperature. 


Since  gas  temperature  in  both  core  flow  and,  with  proper  active  wall  cooling,  near  electrodes  is  below  1000  K, 
thermal  ionization  is  also  negligible.  Therefore,  electron  beam-sustained  plasmas  in  hypersonic  MHD  devices 
should  have  higher  thresholds  of  ionization  and  thermal  instabilities  than  conventional  thermally  ionized  plasma 
would.  We  emphasize  that  the  intersegment  arcing  in  cold-air  MHD  devices  with  electron  beam-controlled 
ionization  was  never  studied,  and  it  should  be  studied  in  the  future  because  of  the  critical  importance  of  E^  for  the 


performance  of  those  devices. 

The  analysis  in  the  preceding  paragraphs  of  this  subsection  was  so  far  focused  on  MHD  interaction 
parameters  in  the  core  flow.  Similarly  to  the  EHD  analysis  in  subsection  2.2,  one  can  define  an  MHD  interaction 
parameter  with  respect  to  shear  stress  at  the  wall; 


aB^L 


pu 


(25) 


where  u'  is  the  friction  velocity  defined  by  Eq.  (8).  If,  for  rough  estimates,  we  assume  that  the  conductivity  cr  and 
the  density  p  near  the  wall  are  equal  to  those  in  the  core  flow  (this  can  be  ensured  by  a  proper  wall  cooling  and  by 
contouring  the  profiles  of  ionizing  electron  beams),  then  from  Eqs.  (25),  (8)-(10)  we  obtain: 


^  (jB^L 

S^=  — X - 

.  ■y  pu 


(26) 


Therefore,  in  the  range  Re’  =10* -10’,  the  MHD  interaction  parameter  with  respect  to  wall  shear  stress  greatly 
exceeds  5:  x  305  .  This  may  have  interesting  implications  for  wall  friction  and  transition  control. 

Of  course,  MHD  effects  on  the  boundary  layer  would  be  quite  complex.  First,  if  the  load  factor  k  is  not 
close  to  1  or  0,  then  the  Joule  heating  is  characterized  by  the  same  interaction  parameters  (15)  and  (25),  and  the 

heating  effects  should  be  comparable  to  those  of  ponderomotive  jy.B  forces  both  in  the  core  flow  and  near  the 
wall.  Second,  if  the  MHD  region  extends  into  the  core  flow,  the  changes  in  the  core  flow  would  also  affect  the 
boundary  layer.  Thus,  evaluation  of  MHD  effects  on  the  boundary  layer  requires  a  fully  coupled  analysis. 


III.  Electrodeless  MHD  flow  control:  a  magneto-thermal  funnel 


36 


Since  electrode-related  problems  like  intersegment  arcing  can  potentially  reduce  performance  of  MHD 
flow  control  devices,  it  is  interesting  to  consider  whether  MHD  flow  turning  or  compression  can  be  accomplished 
without  electrodes. 

Using  one  of  the  Maxwell’s  equations 

VxB  =  MoJ  (27) 


and  the  vector  identity 


we  obtain  from  Eq.  (5): 

p—  =  -Vp  +  e{n^ -n_  ~nAE  +  —(B-V)B-'V^—+pAu 
At  '  n  \ 


Bo 


=  -V 


P  + 


2m 


2  ^ 


+  — +  -n^')E+fuAu 


0/ 


(28) 


(29) 


8^  1 

Thus,  magnetic  field  energy  density,  - ,  acts  as  an  additional  pressure.  The  term  — (B-ViB  can  be 

2/^0  Bo 

interpreted  as  tension  of  magnetic  field  lines  due  to  their  curvature.  These  ponderomotive  forces  are  successfully 
used  for  containment  of  fully  ionized  plasmas  in  fusion  devices. 

With  S=1  Tesla,  the  magnetic  pressure  is  about  8  atm,  and  significant  MHD  flow  turning  and  compression 
effects  might  be  expected  in  hypersonic  flow.  However,  practical  applications  of  the  magnetic  ponderomotive  flow 
control  would  be  hampered  by  low  ionization  level.  Indeed,  with  the  conductivity  of  only  0.1-1  mho/m  typical  for 
MHD  devices  with  ionization  by  electron  beams,  the  magnetic  Reynolds  number  Re^  =  p^auL  □  1 ,  and  the 
magnetic  field  is  not  “frozen”  into  the  plasma.  Therefore,  ponderomotive  forces  depend  on  the  conductivity,  and 
their  relative  strength  is  characterized  by  the  Stuart  number  5  that,  as  shown  in  the  previous  section,  is  also  less  than 
1  for  an  MHD  region  of  a  few  meters  in  length.  Note  the  important  difference  between  electric  and  magnetic 
pressures  in  weakly  ionized  plasmas.  While  the  electric  pressure  can  fully  act  on  even  weakly  ionized  gases  (see 
section  2.2),  the  magnetic  pressure  cannot  be  fully  exerted  on  weakly  ionized,  low-conductivity  media,  unless  the 
flow  velocity  and/or  the  length  scale  are  extremely  large. 

Consider  a  sample  case  of  electrodeless  MHD  flow  control,  illustrated  in  Fig.  2.  Gas  flow  into  a  cylindrical 
duct  (inlet),  and  magnetic  field  lines  are  coming  out  of  the  duct  in  a  pattern  similar  to  that  near  the  edge  of  a 
solenoid.  Ionizing  electron  beams  form  an  annular  ring  coaxial  with  and  adjacent  to  the  duct  wall,  and  are  injected 
upstream,  along  B  field  lines.  With  the  plasma  moving  across  the  nonuniform  magnetic  field,  induced  electric 

currents  form  circles  coaxial  with  the  inlet.  The  resulting  jxB  forces  decelerate  and  compress  the  flow. 
Additionally,  gas  heating  by  the  induced  currents  creates  reduced-density  region  around  and  upstream  of  the  inlet, 
deflecting  streamlines  and  causing  an  increase  in  compression  and  mass  flow  rate  into  the  duct.  Thus,  both  jxB 
forces  and  heating  create  what  may  be  called  a  magneto-thermal  funnel. 

The  calculations  were  performed  for  ffeestream  conditions  corresponding  to  Mach  8  flight  at  the  altitude 
A=30  km,  with  static  pressme  and  temperature  of  1 197  Pa  and  226.5  K,  respectively.  The  cylindrical  inlet  of  radius 
'into  =  1  in  was  positioned  at  axial  location  =  2  m.  The  model  profile  of  a  solenoid-like  stationary  magnetic 
field  in  the  upstream  region  was: 

B=[B,(r,z),B,(r,z),0];  B,(r,z)  =  -(s’(z,r)-R,^(r,z)f'; 

=  5a.n,»  exp(-(z  -  r,.„^)/4)  ;  5, (r, z)  =  R^(z)[l  -  exp(-(r  /(0.5  •  ?;„,,,))] ,  (30) 


with  =  7  T. 

The  ionizing  electron  beams  are  injected  in  a  1-cm  wide  ring  around  the  inlet.  The  energy  of  beam 
electrons  is  ®  50  keV,  and  their  relaxation  length  under  the  conditions  of  this  case  is  «  2.5  m.  The  current 
density  of  the  beam  at  the  injection  point  is  very  high,  «  72  mA/cm^,  and  the  total  power  of  the  beam  is  Qj=4.57 
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MW.  Power  deposition  by  electron  beams  was  modeled  by  a  2D  profile  approximately  following  magnetic  field 
lines: 

9, (z,r)  =3x10“' expH(z-0.7)V0.45^+(r-rjVr^(zf]}  W/m\ 

=1.75exp(-z/0.5)  +  0.967;  r;^(z)  =  0.01 +  0.2 exp(-z/ 1.05) .  (31) 

where  z  and  r  are  in  meters.  The  profile  of  both  beam-induced  ionization  rate  corresponds  to  the  beam  power 
deposition  profile: 

(32) 

where  W.  =  34  eV  is  the  energy  cost  of  ionization  by  high-energy  beam  electrons. 


The  generalized  Ohm’s  law  (13)  in  the  case  u  =  («.,«,, 0) ,  B  =  ,  j  =  (0,0,  j^)  is: 

crEl  _a(u^B^-B^u^) 


I  = - ^ 

’’  l  +  Op,  1  +  Qp,. 

This  azimuthal  current  produces  Joule  heating  of  the  gas: 

and  jxB  forces: 

F.=-  jB^  and  F  =  jB. . 

z  J  q)  r  T  J  (p  z 

The  fiill  set  of  Euler  equations  in  cylindrical  coordinates  is: 

—U  +  —  R^—Z  =  H, 
dt  dr  dz 
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p  =  (y-l)p5j„, 
e,o,=F[^in,+(«/+“/)/2], 


(33) 

(34) 

(35) 

(36) 
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(38) 

(39) 


where  p  and  p  denote  the  gas  density  and  pressure,  respectively,  are  the  r  and  z  velocity  components,  e,oi  is 
the  total  energy  (internal  energy,  excluding  the  energy  of  vibrational  mode,  plus  kinetic  energy)  of  the  gas  per  unit 
volume,  and  is  the  internal  energy  (excluding  the  energy  of  vibrational  mode)  per  unit  mass. 

The  source  term  in  the  total  energy  equation  is: 

Q  =  -Q„  +  Qvt  +  >  (40a) 

where  the  vibrational  excitation  term  can  be  expressed  as  a  fraction  r\^  of  the  Joule  heating  rate:  =  77^7^  la . 

The  fraction  rj^  is  a  function  of  the  local  reduced  electric  field  =  \u^B,.-B^u)^ln .  The  fraction  rj^  was  taken 

from  Ref.  49,  where  it  was  tabulated  as  a  function  of  E In  on  the  basis  of  solution  of  Boltzmann  kinetic  equation 
for  plasma  electrons.  The  remaining  energy  addition  and  dissipation  terms  are: 

Qy,={E^-EyT))lT„(T),  (40b) 
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heating  rate  per  unit  volume  due  to  VT-relaxation;  and  -  power  deposited  per  unit  volume  by  the  electron  beam. 
Nonequilibrium  and  equilibrium  vibrational  energy  was  expressed  through  the  respective  temperatures  by  the  Planck 
formula,  and  nitrogen  vibrational  relaxation  time  was  taken  as  in  Refs.  16-18. 

The  plasma  was  modeled  as  consisting  of  electrons,  positive  and  negative  ions,  whose  number  densities 
obey  the  quasineutrality:  »  n^  +«_ .  The  set  of  equations  for  kinetics  of  charge  species,  accounting  for 
electron  beam-induced  ionization  rate  term),  ionization  rate  due  to  plasma  electrons  (with  Townsend  ionization 
coefficient  a  depending  on  the  local  reduced  electric  field,  - Bju^jn"),  attachment  of  electrons  to 

molecules  with  formation  of  negative  ions  (frequency  v^),  collisional  detachment  of  electrons  from  negative  ions 
(rate  constant  kj),  and  electron-ion  and  ion-ion  recombination  (rate  coefficients  and  fin,  respectively),  jg;i«-'*'33,36 

-^+divf^=a\j^\/e+q,+kjNn_ -vji^ 

-^^divf^=a\jj^le  +  q^-p„n_n^  (41) 

^  +  divf_  =  -kjNn_  + 
at 

In  the  considered  case,  there  are  no  drift  fluxes  along  either  z  or  r  axes  (only  azimuthal  current,  ,  is  present),  and. 


therefore,  the  fluxes  of  charged  species  in  (r,  z)  plane  can  be  written  simply  as  E ^ ^  _(r,z)  =  +u^  -i^) , 


where  i)  .  are  the  vectors  of  unit  length  in  z  and  r  directions; 
magnetic  field  taking  into  accoimt  ion  slip,”’ 


i+Q,n,  ■ 


is  the  effective  electron-ion  velocity  across 


(42) 


The  initial  conditions  for  plasma  components  are:  +  .(r,z,t  =  0)  =  0  .  The  boundary  conditions  are: 
n,..._(z  =  0,r,/)  =  (z  =  z„„,r,t)  =  0 

(43) 

=  r^,t)  =  0;  dn^^_{z,r  =  0,/)/dr  =  0. 

Rate  coefficients  of  electron-ion  and  ion-ion  recombination,  and  of  electron  attachment  and  detachment 
processes,  discussed  in  Refs,  16-18,  31-34,  were  taken  from  Refs.  50  and  51.  Since  some  of  those  rate  coefficients 
depend  on  electron  temperature  T^,  it  is  important  to  calculate  the  electron  temperature  in  the  modeling.  In  our 
computations,  electron  temperature  was  determined  from  the  tabulated  data  on  electron  diffusion  and  mobility 
coefficients  of  Ref.  49.  In  that  paper,  the  diffusion  and  mobility  coefficients  are  listed  as  functions  of  E/N, 
determined  from  experimental  data  and  extrapolation  based  upon  solution  of  Boltzmaim  kinetic  equation  for 
electrons  in  air. 

Results  of  the  calculations  are  shown  in  Figs.  3-10.  Fig.  3  depicts  the  magnetic  field  lines,  and  Fig.  4  shows 
the  electron  beam  power  deposition  profile.  The  resulting  profiles  of  electron  density  and  azimuthal  current  are 
shown  in  Figs.  5  and  6.  Figs.  7-10  show  the  computed  profiles  of  static  and  vibrational  temperatures,  gas  density  and 
radial  velocity.  The  formation  of  magneto-thermal  funnel  is  clearly  seen  in  the  figures. 

The  overall  performance  of  the  fuimel  can  be  assessed  by  the  mass  flow  rate  and  enthalpy  flux  into  the 
inlet,  and  by  the  total  pressure  coefficient.  Without  ionization  ( ^^  =  0  ),  the  mass  flow  and  enthalpy  flux  into  the 


inlet  are  the  following:  m_.  ^=  ^2nrp{z,r)u.dr  =139.6  kg/s,  and 

0 


^  =  J 2n:rp(z,r)u.lyp/(}' -\)p  +  0.5{ul  +u^)]dr  =438.5  MW. 

0 

The  total  pressure  coefficient,  that  is,  the  total  pressure  averaged  over  the  inlet  cross-section  divided  by  the 
freestream  total  pressure,  is  <  >1 P, a o~^-  When  the  electron  beam  is  turned  on,  but  in  the  absence  of  any 

magnetic  field  (B=0,  g ^  ^  0 ),  gas  heating  by  the  beam  and  streamline  deflection  result  in  a  barely  noticeable 
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increase  in  mass  and  enthalpy  fluxes,  and  a  slight  decrease  in  total  pressure:  =  139.9  kg/s;  ^  =439.4 

MW;  <  j  =0.997.  With  both  ionization  and  magnetic  field  on  ( .5  =  B{r,z) ,  the  increased  mass 

and  enthalpy  fluxes  and  the  decreased  total  pressure  coefficient  are:  =146  kg/s;  ^  =457.8  MW; 

<  Plot  > ! Piot.o  =0-93.  The  maximum  flow  turning  induced  in  the  flowfield  is  about  2.5°. 

The  computations  were  performed  with  2“'*  order  MacCormack  method^^  on  a  rectangular  310x250  grid. 
Computations  with  a  finer  grid,  465x375,  show  that  the  mass  and  enthalpy  fluxes  and  the  total  pressure  are  within 
0.1%  from  those  computed  with  310x250  grid,  thus  confirming  the  accuracy  of  computations. 

Thus,  magneto-thermal  funnel  can  indeed  increase  mass  capture  and  total  enthalpy  flux  into  the  inlet. 
However,  because  of  gas  heating,  the  total  pressure  would  decrease.  In  the  sample  case,  the  beam  and  magnetic  field 
parameters  were  not  optimized,  and  the  effects  are  quite  modest,  despite  the  very  high  current  density  of  the  beam. 
The  mass  flow  rate  increase  can  be  made-  stronger  by  using  stronger  magnetic  fields  and  higher  beam  currents. 
Practicality  of  this  method  of  flow  control,  with  advantages  weighted  against  flaws  (strong  magnetic  fields,  high 
electron  beam  currents,  and  losses  of  total  pressure)  should  be  a  subject  of  a  systems  study. 

IV.  Conclusions 

The  principal  difficulty  in  high-speed  flow  control  using  electric  and  magnetic  fields  is  that  the  relatively 
cold  gas  has  to  be  ionized  in  electric  discharges  or  by  electron  beams,  which  requires  large  power  inputs  and  results 
in  low  ionization  fraction  and  electrical  conductivity.  The  low  ionization  fraction  means  that,  although  electrons  and 
ions  can  interact  with  electromagnetic  fields,  transfer  of  momentum  and  energy  to  or  from  the  bulk  neutral  gas  can 
be  small  compared  with  momentum  and  energy  carried  by  the  high-speed  flow. 

Even  at  the  highest  values  of  the  electric  field  that  can  exist  in  cathode  sheaths  of  electric  discharges, 
electrohydrodynamic  (EHD),  or  ion  wind,  effects  in  a  single  discharge  can  be  of  significance  only  in  low-speed  core 
flows  or  in  laminar  sublayer  of  the  boundary  layer.  Cumulative  action  of  multi-element  discharges  would 
conceivably  amplify  the  single-sheath  effect.  However,  Joule  heating  can  overshadow  the  cathode  sheath  ion  wind 
effects. 

Theoretical  analysis  of  MHD  flow  control  with  electron  beam  ionization  of  hypersonic  flow  shows  that  the 
MHD  interaction  parameter  is  a  steeply  increasing  function  of  magnetic  field  strength  and  the  flow  velocity. 
However,  constraints  imposed  by  arcing  between  electrode  segments  can  reduce  the  performance  and  make  the 
maximum  interaction  parameter  virtually  independent  of  Mach  number.  Thus,  the  value  of  critical  Hall  field 
corresponding  to  the  onset  of  arcing  between  electrode  segments  is  extremely  important  in  designing  and  evaluating 
performance  of  hypersonic  MHD  devices,  and  the  arcing  instability  should  be  investigated  theoretically  and 
experimentally. 

Estimates  of  MHD  interaction  parameter  rvith  respect  to  shear  stress  show  that  the  relative  strength  of 
MHD  effects  can  be  higher  near  the  wall  (in  the  boundary  layer)  than  in  the  core  flow,  which  may  have  implications 
for  MHD  boundary  layer  and  transition  control.  Overall,  one  of  the  qualitative  conclusions  from  this  work  is  that 
large-scale  bulk  MHD  effects  on  the  overall  hypersonic  flow  field  are  quite  difficult  (although  not  impossible)  to 
achieve,  whereas  small  and/or  localized  MHD  effects  such  as  those  for  small-angle  flow  turning  or  boundary  layer 
control  may  be  useful  for  significant  flow  field  manipulation. 

In  this  work,  we  also  suggested  a  magneto-thermal  funnel  concept  that  would  accomplish  and  electrodeless 
MHD  turning  and  compression  of  high-speed  flows.  Computations  of  a  sample  case  demonstrate  that  the  turning  and 
compression  of  hypersonic  flow  ionized  by  electron  beams  can  be  achieved;  however,  the  effect  is  relatively  modest 
due  to  low  ionization  level. 
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Fig.  3.  Magnetic  field  lines  in  the  computed  case. 
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Fig.  4.  Electron  beam  power  deposition  contour  lines  in  the  computed  case 


Fig.  5.  Electron  number  density  contours  in  the  computed  case. 
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Fig.  7.  Static  temperature  contours  in  the  computed  case. 


Fig.  8.  Vibrational  temperature  contour  lines  in  the  computed  case 
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Fig.  10.  Radial  velocity  contour  lines  in  the  computed  case 
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CHAPTER  3 

NONEQUILIBRHJM  MHD  CONTROL  OF  SCRAM  JET  INLETS 


Summary 

This  chapter  is  devoted  to  theoretical  analysis  of  magnetohydrodynamic  (MHD)  control  of  forebody  flow 
compression  and  air  mass  capture  in  scramjet  inlets  for  vehicles  that  would  fly  at  Mach  6-10.  Due  to  the  low  static 
temperature,  nonequilibrium  air  ionization  and  electrical  conductivity  are  created  by  electron  beams  injected  into  the 
gas  along  magnetic  field  lines.  Two-dimensional  inviscid  steady-state  flow  equations  are  solved  jointly  with 
equations  describing  electron  beam-induced  ionization  profiles,  plasma  kinetics,  and  MHD  equations.  Several 
scenarios  are  considered.  At  flight  Mach  numbers  8  and  10,  with  forebody  and  inlet  geometry  designed  for  Mach  6, 
the  shocks  that  would  enter  the  inlet  can  be  moved  back  on  the  cowl  lip  by  placing  an  MHD  generator  at  the  first 
compression  ramp.  In  addition  to  controlling  the  shocks,  the  MHD  device  would  generate  electric  power;  however. 
Joule  heating  would  result  in  losses  of  total  pressure.  Another  scenario  is  an  on-ramp  MHD  accelerator  that  should 
increase  mass  capture  at  Mach  6  or  8  for  a  vehicle  designed  for  Mach  1 0.  Calculations  show  that  this  scenario  has 
only  disadvantages,  as  the  MHD  device  consumes  high  power,  reduces  total  pressure,  and  actually  decreases  mass 
capture  due  to  Joule  heating  and  thermal  expansion  of  the  gas.  A  modest  increase  in  mass  capture  can  be  in  principle 
accomplished  in  an  MHD  generator  mode,  if  the  magnetic  field  has  components  both  parallel  and  orthogonal  to  the 
flow.  However,  this  scenario  requires  unrealistically  large  volumes  of  strong  magnetic  fields,  and  the  mass  capture 
increase  is  due  mostly  to  a  nonuniform  gas  heating.  A  new  concept  of  a  virtual  inlet  is  proposed,  where  a  localized 
off-body  energy  addition  is  used  to  increase  mass  capture,  while  not  reducing  (and  even  increasing)  total  pressure  at 
the  inlet. 

1.  Introduction 

The  geometry,  size,  and  weight  of  scramjet-powered  hypersonic  vehicles  are  largely  dictated  by  the  need  to 
compress  the  ambient  low-density  air  upstream  of  the  combustor.  The  optimum  geometry  corresponds  to  the  well- 
known  shock-on-lip  (SOL)  condition  (Fig.  1):  the  compression  ramp  shocks  converge  on  the  cowl  lip,  and  the 
reflected  shock  impinges  on  the  upper  boundary  of  the  inlet.  Since  shock  angles  are  determined  by  the  flight  Mach 
number,  the  SOL  condition  cannot  be  met  at  Mach  numbers  higher  or  lower  than  the  design  Mach  number.  At  Mach 
numbers  higher  than  the  design  one,  the  shocks  move  inside  the  inlet  (Fig.  2),  causing  multiple  reflected  shocks,  loss 
of  stagnation  pressure,  possible  boundary  layer  separation,  and  engine  unstart.  At  Mach  numbers  lower  than  the 
design  one,  the  so-called  spillage  occurs,  and  the  air  mass  capture  decreases  (Fig.  3). 

To  avoid  performance  penalties  at  off-design  Mach  numbers,  a  variable  geometry  inlet  can  be  used.  However,  the 
mechanical  variable  geometry  system  would  be  quite  heavy,  weighing  perhaps  tens  of  metric  tons.  An  alternative 
approach  is  to  optimize  inlets  using  energy  addition  to  or  extraction  from  the  flow.  Plasmas  and  various 
magnetohydrodynamic  devices  may  offer  a  viable  optimization  scheme. 

One  method  that  was  analyzed  by  several  groups,'  *  is  to  design  the  vehicle  for  a  relatively  low  Mach  number  (for 
example,  Mach  6),  and  move  the  shocks  from  inside  the  inlet  back  to  the  cowl  lip  by  an  MHD  generator  device 
placed  at  one  of  the  compression  ramps  (Fig.  4).  The  advantages  of  this  method  are  the  following:  1)  reduction  in 
vehicle  size  and  weight  due  to  both  absence  of  mechanically  variable  inlet  and  the  shortening  of  compression  ramps 
because  of  the  reduction  of  design  Mach  number;  the  weight  savings  may  more  than  offset  the  weight  of  the  magnet 
and  other  MHD-related  hardware;  2)  no  net  power  is  required  to  mn  the  MHD  device,  since  the  generator  mode  is 
used,  and  the  power  requirements  for  air  ionization  can  be  minimized  with  electron  beams  as  ionizers;  3)  flexibility 
of  flow  control.  An  important  disadvantage  of  this  approach  is  that  the  inevitable  Joule  heating  and  other  dissipative 
processes  result  in  substantial  stagnation  pressure  losses. 

Another  method  of  MHD  control  is  to  design  the  vehicle  for  a  high  Mach  number  (for  example,  Mach  10),  and  to 
decrease  or  eliminate  the  flow  spillage  by  an  MHD  accelerator  device  placed  at  one  of  the  compression  ramps  (Fig. 
5).  A  substantial  disadvantage  of  this  method  compared  with  the  MHD  generator  method  shown  in  Fig.  4  is  that  an 
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accelerator  requires  a  substantial  amount  of  power.  Meanwhile,  the  accelerator  method  would  still  result  in 
considerable  stagnation  pressure  losses. 

Yet  another,  recently  suggested,^  method  of  MHD  inlet  control  is  to  operate  an  on-ramp  MHD  device  in  generator 
mode,  but  use  a  magnetic  field  configuration  with  B  field  components  both  parallel  and  orthogonal  to  the  flow. 
According  to  Ref.  5,  air  mass  capture  at  Mach  numbers  lower  than  the  design  one  can  be  increased  with  this  method. 

In  this  paper,  we  perform  computational  and  qualitative  analysis  of  all  the  three  methods  of  MHD  flow  and  shock 
control  at  off-design  Mach  numbers.  We  also  propose  a  new  method  of  inlet  control:  a  virtual  inlet. 

An  important  factor  in  the  analysis  is  that  in  the  flight  regime  of  interest,  static  temperature  is  too  low  for  thermal 
ionization  of  air.  Therefore,  nonequilibrium  ionization  in  MHD  devices  must  be  done.  The  energy  cost  of  ionization 
imposes  very  rigid  constraints  on  the  choice  of  ionization  methods:  only  the  most  energetically  efficient  ionization 
method  can  be  used  in  these  devices.  Our  earlier  work  showed  that  electron  beams  injected  into  the  gas  along 
magnetic  field  lines  constitute  the  most  efficient  way  of  creating  nonequilibrium  ionization.  However,  even  with  the 
most  efficient  ionizer,  the  very  existence  of  ionization  cost  does  not  allow  the  ionization  degree  and  the  conductivity 
in  hypersonic  MHD  devices  to  be  high.  This  makes  the  task  of  using  MHD  for  hypersonic  flow  control  and  power 
generation  quite  challenging. 

II.  The  model 

We  consider  hypersonic  gas  flow  along  the  two  ramps  with  10®  and  20®  angles  upstream  of  the  inlet  with  forward- 
shifted  cowl  lip,  as  shown  in  Figs.  1-6.  The  flow  is  two-dimensional  in  (x,  z)  plane.  Cases  both  without  and  with 

MHD  influence  on  the  flow  are  computed.  In  MHD  cases,  both  magnetic  field  and  ionizing  electron  beam  are 
directed  parallel  to  z-axis.  Because  the  entire  flow  region  is  hypersonic,  steady  state  solution  using  x  as  marching 
coordinate  can  be  found. 

The  set  of  Euler  equations  in  Cartesian  coordinates,  together  with  an  ideal  gas  equation  of  state,  a  simple  model  of 
an  ideal  Faraday  MHD  generator,  vibrational  relaxation,  and  a  plasma  kinetic  model  are  those  of  Refs.  7  and  8.  In 
the  present  paper,  we  solve  steady-state  problems  only,  so  that  all  time  derivatives  in  the  equations  of  Ref.  8  were 
equated  to  zero.  Additionally,  in  MHD  accelerator  cases,  the  load  factor  k>\,  so  that  in  the  equations  of  Ref.  8, 
(1-k)  should  be  substituted  by  -1) ,  and  the  power  P  dissipated  in  the  external  load  per  unit  gas  volume  should 
be  taken  with  the  sign  opposite  to  its  sign  in  Eq.  (5)  of  Ref.  8.  In  the  cases  with  B  field  with  both  B^  and 

components,  the  corresponding  components  of  y  x  R  force,  as  well  as  both  contributions  to  the  Faraday  e.m.f.,  were 
included  in  the  equations,  similar  to  Ref.  5. 

Power  deposition  and  ionization  by  the  electron  beams  have  to  be  coupled  with  gas  dynamic  and  plasma  kinetic 
equations.  In  our  earlier  work,®  '®  we  used  the  so-called  “forward-backward”  approximation'®  for  electron  beam 
propagation  along  magnetic  field  lines  into  gases  and  beam-generated  ionization.  This  method  is  quite  accurate,  and 
it  can  be  (and,  in  fact,  was®’  *°)  coupled  with  Navier-Stokes  and  MHD  equations.  However,  computational  demands 
of  the  fiill  “forward-back”  method,  when  coupled  with  gas  dynamic  and  MHD  equations,  are  prohibitive.  Therefore, 
a  simple  but  accurate  way  to  predict  ionization  profiles  should  be  found. 

Having  analyzed  physics  of  ionization  processes,  we  concluded  that  in  a  uniform  gas  the  ionization  rate  profile  has 
to  be  close  to  a  tmncated  Gaussian.  Indeed,  ionization  of  molecules  by  electron  impact  is  inefficient  both  at  low  and 
high  (above  1  keV)  electron  energies,'®’ '®  The  maximum  of  the  ionization  cross  section  corresponds  to  electron 
energy  of  several  hundred  electronvolts.'®’  In  this  energy  range,  and  even  above  it,  the  probability  of  back 
scattering  is  substantial,  and  there  are  both  forward  and  backward  electron  fluxes  exchanging  electrons  between 
them.  The  resulting  motion  is  very  much  like  diffusion.  Thus,  energy  deposition  (Q,)  and  ionization  profiles 
computed  with  both  “forward-back”  method  and  CYLTRAN  Monte  Carlo  Calculations  (MCC)®''®  are  very  close  to 
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Gaussian  profiles:’’®  =  a  +—exp(-2(^  ^ 

w 

(1) 

where  <^  =  Zi^(x)-z;  z^«I^/3.21,  wail.64z„;  L^(e^,N)  is  the  beam  relaxation  length.  There  are  two  additional 
conditions  for  determining  the  constants  a  and  b: 

Q,{L,)  =  0  and  (2) 

0 

£•(,  and  Ji^  are  the  beam  electron  energy  and  current  density  at  the  injection  point. 

The  use  of  a  loss  function’’®  or  the  direct  solution  of  kinetic  equation  for  electrons  in  the  “forward-back” 
approximation  give  a  relaxation  length  overestimated  by  a  factor  of  1.5-2.  This  is  because  lateral  electron  scattering 
and  the  loss  of  an  additional  energy  equal  to  the  energy  of  the  secondary  electron  were  neglected  in  those 
approximations.  In  the  present  paper,  an  empirical  approximation  formula  for  the  electron  beam  relaxation  length 
that  effectively  includes  those  two  effects  is  used. 

Assuming  that  the  empirical  approximation*’  for  the  beam  relaxation  length  in  a  gas  with  uniform  number  density, 
A=const,  is  also  valid  for  nonuniform  neutral  gas  density,  but  with  the  average  number  density  along  the  beam 
direction,  N ,  and  that  the  relaxation  length  along  magnetic  field  in  dense  gases  is  the  same  as  without  the  magnetic 
field,  we  have  for  the  relaxation  length: 

4=l.M0”fi'’//7,m  (3) 

where  is  in  keV,  N  is  in  m'®,  and  Z,^  is  in  meters.  At  each  location  along  x,  with  the  Ij,  in  this  location,  we  can 

find  N(x)  and  : 

1  ^/t 

N(x)  =  -^  jNi^,x)d4;  £,(x)  =  (I,A/l.MO’*)''”,keV.  (4) 

0 

Equations  (l)-(4)  fully  determine  both  beam-induced  profiles:  energy  deposition,  and  ionization  rate: 
qAx,4)«Q,{x,^)l{eW.:),  (5) 

where  W.  =  34  eV  is  the  energy  cost  of  ionization  by  high-energy  beam. 

In  Fig.  7,  obtained  with  Eqs.  (l)-(4)  is  compared  with  results  of  CYLTRAN  Monte  Carlo  calculations  (MCC).^’*® 
Note  that  when  the  Gaussian  approximation  (1)  uses  the  same  relaxation  length  as  that  given  by  the  MCC,  the 
agreement  is  very  good. 

III.  Computed  cases 

MHD  control  at  M>M(iesjg„ 

The  inlet  geometry  was  first  optimized  for  Mach  6:  the  location  of  the  cowl  lip  was  chosen  so  that  both  oblique 
shocks  would  together  reach  the  lip,  and  the  reflected  shock  would  fall  right  on  the  comer  between  the  second  ramp 
and  the  engine.  This  configuration  is  shown  in  Fig.  8.  The  ffeestream  conditions  correspond  to  the  altitude  of  26  km. 
The  shocks  inside  the  inlet  are  weak,  as  can  be  seen  from  the  streamline  pattern. 

Next,  the  flow  velocity  was  increased  to  Mach  8  for  the  same  geometry.  To  keep  the  freestream  djuiamic  pressure 
about  the  same  as  at  Mach  6,  the  flight  altitude  has  to  be  increased,  and  all  Mach  8  cases  were  ran  for  an  altitude  of 
30  km,  and  Mach  10  cases  -  for  33  km.  As  seen  in  Fig.  9,  the  shock  would  now  make  contact  ■with  the  cowl  farther 
downstream,  producing  strong  multiple  reflected  shocks.  Although  in  this  work,  shock-boundary  layer  interactions 
and  unsteady  processes  were  not  studied,  flow  separation  and  unstart  can  be  expected  in  this  case. 

To  find  out  if  an  MHD  system  could  bring  the  flow  closer  to  the  design  point,  an  ideal  Faraday  generator  system 
with  the  loading  parameter  k  =  0.5  was  put  near  the  surface  of  the  first  ramp.  Both  the  magnetic  field  and  the 
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electron  beam  were  directed  downward  along  z-axis.  Both  the  B  field  and  the  beam  were  assumed  uniform  along  x 
within  the  region  2.75  m<x<  10,8  m  ,  and  to  be  equal  to  zero  elsewhere.  The  presumed  MHD  configuration  in 
(x,  z)  and  (y,  z)  planes  is  shown  schematically  in  Fig.  6.  The  B  field  profile  along  the  z-axis  was  Gaussian 


(0.75  m)" 


with  the  maximum  of  B^  =  4  Tesla  at  the  wall.  Electron  beam  energy  and  current 


K  '  '  j 

density  were  adjusted  so  as  to  move  the  shocks  at  Mach  8  back  to  the  cowl  lip.  In  this  case  and  throughout  this 
paper,  the  electron  beam  current  was  set  at  =  1  mA/cm^  . 


As  seen  in  Fig.  10,  with  the  MHD  generator  on,  both  shocks  can  be  moved  back  to  the  lip,  and  streamlines  in  the 
inlet  duct  become  much  more  parallel  to  the  walls.  However,  their  angles  of  incidence  are  slightly  different  from 
those  in  the  Mach  6  case,  so  that  the  reflected  shock  would  not  exactly  follow  the  design  line.  Perhaps  an  additional 
adjustment  with  an  internal  MHD  or  other  control  system  would  be  required. 

Details  of  charged  species  and  vibrational  temperature  profiles  in  the  computed  MHD  case  are  shown  in  Fig.  1 1.  As 
seen  in  the  figure,  inside  the  MHD  region  with  electron  beam  ionization,  electron  density  is  of  the  order  of 
10'^  cm‘^ ,  controlled  by  beam-induced  ionization  and  electron-ion  recombination.  Attachment  is  insignificant,  and 
the  negative  ion  density  is  small.  Downstream  of  the  beam-ionized  MHD  region,  electrons  disappear  almost 
immediately,  but  the  wake  of  about  5x10*  cm  *  negative  ions  persists  and  enters  the  inlet.  Vibrational  temperature 
(Fig.  11)  does  exceed  the  gas  temperature  (Fig.  10),  but  the  degree  of  nonequilibrium  and  the  enthalpy  content  of  the 
vibrational  mode  are  not  high. 

At  Mach  10,  h=33  km  (design  regime  M=6,  h=26  km)  (Fig.  12a,  b),  return  of  the  shocks  back  to  the  cowl  hp 
requires,  with  the  same  electron  beam  current  density  and  energy  and  the  same  maximum  B  field,  an  adjustment  to 
the  length  and  location  of  the  MHD  region  compared  with  fliat  in  the  Mach  8  case.  The  MHD  region  is  now 
positioned  at  4  m<x<  12.25  m  .  Note  that  in  principle  a  change  of  MHD  region  location  and  size  can  be  achieved 
by  controlling  the  electron  beam  profile,  without  moving  or  changing  the  magnets.  In  Fig.  13,  electron  beam  energy 
at  the  injection  point  is  plotted  versus  distance  along  the  MHD  generator  region  in  the  cases  of  Fig.  10,  11,  and  12, 
with  =2.3  m  =  const. 

Table  1  lists  some  important  parameters  in  the  computed  cases  at  the  inlet  entrance  cross  section  (at  the  location 
corresponding  to  the  coiner  at  the  end  of  the  last  ramp):  mass  ( m)  and  total  enthalpy  flow  rate  {H)  into  the  inlet, 
and  averaged  over  the  inlet  cross  section  compression  ratio,  total  pressure  recovery  coefficient,  translational 
temperature,  and  mass  capture.  As  seen  in  the  Table,  MHD  operation  is  accompanied  by  substantial  losses  of  total 
pressure.  Table  2  lists  MHD  generator  parameters:  extracted  power,  ,  electron  beam  input  power,  ,  power  of 

the  y  X  5  force,  P-.^^ ,  Joule  dissipation  rate,  Pj ,  the  difference  between  the  MHD-generated  electric  power  and  the 
beam  power,  (P^Hp  —Pi,),  and  enthalpy  extraction  ratio  (i.e.,  the  MHD-generated  electric  power  as  a  percentage  of 
the  total  enthalpy  entering  the  inlet  per  unit  time),  x  ■ 


The  enthalpy  extraction  ratio  in  both  MHD  cases  is  about  2-4%.  The  power  cost  of  ionization  is  quite  low  compared 
with  other  power  budget  components,  including  the  generated  electric  power.  Therefore,  the  ionizing  beams  can  be 
powered  by  a  fraction  of  generated  power,  and  a  substantial  portion  of  that  power  can  be  used  for  other  purposes. 
The  losses  of  total  pressure  seen  in  Table  1  are  due  to  entropy-generating  Joule  dissipation. 

All  the  computed  cases  were  inviscid.  We  have  performed  a  preliminary  uncoupled  estimate  of  the  boundary  layers 
for  design  Mach  6,  h=26  km  and  off-design  Mach  8,  h=30  km  cases  with  and  without  the  MHD  at  x=10  m 
downstream  from  the  nose.  The  boundary  layer  was  modeled  as  that  for  viscous  flow  along  the  flat  plate,  using  the 
steady  state  core  flow  profiles  of  pressure,  velocity,  temperature  and  current  density  (in  the  MHD  case)  computed 
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with  the  inviscid  model  as  boundary  conditions.  In  the  MHD  case,  the  jxB  body  force  acting  on  the  flow  in  the 
upstream  direction,  /  =  -|y^(x)5^|cos0,  was  included  in  the  momentum  equation,  and  power  extraction, 
-\j y(.x)u{x,z)B\^kcosO  -  in  the  energy  equation.  Here  k  =  0.5  is  the  load  factor,  u{x,z)  is  the  flow  velocity  in  the 

boundary  layer,  and  G-\(f  is  the  ramp  angle.  The  wall  temperature  was  assumed  to  be  Twaii=600  K.  The  boundary 
layer  stracture  was  computed  with  of  the  standard  set  of  equations  for  steady-state  compressible  turbulent  flat-plate 
flow'*  developed  by  Cebeci  and  Smith.'®  The  results  are  shown  in  Fig.  14  and  Table  3.  In  Fig.  14  and  Table  3,  5^^ 
is  the  displacement  thickness,  is  the  boundary  layer  thickness  corresponding  to  the  velocity  of  0.99  ffeestream 
velocity,  and  Sj  is  the  thermal  bonndary  layer  thickness  defined  as  the  distance  from  the  wall  where  static 
temperature  is  1.01  of  its  freestream  value. 

b.  MHD  control  at  M<Mdesign:  accelerator  regime 

Consider  the  scenario  when  the  inlet  and  the  forebody  are  designed  for  shock-on-lip  condition  at  a  high  Mach 
number,  and  flow  spillage  occurs  when  flying  at  lower  Mach  numbers  (Fig.  3).  The  task  of  a  plasma  or  MHD  device 
is  to  increase  the  mass  capture. 

The  specific  design  we  consider  is  SOL  at  Mach  10,  h=33  km,  as  shown  in  Fig.  15  (Case  1  in  Table  4).  At  off- 
design  M=6,  h=26  km,  spillage  occurs  (Fig.  16  and  Case  2  in  Table  4).  The  first  method  of  mass  capture  increase  is 
to  use  an  on-ramp  MHD  accelerator  (Fig.  5).  Presumably,  the  accelerator  would  increase  the  compression  ratio  and 
mass  capture.  Figs  .17  and  18  show  results  in  2  accelerator  cases  (Cases  3  and  4  in  Table  4).  In  the  first  one  (Case 
3),  the  magnetic  field  profile  was  identical  to  that  in  the  MHD  generator  case,  and  ionization  by  electron  beams  was 
realistic  with  finite  relaxation  length  L^=2.9  m  and  Gaussian  power  deposition  profile  (Fig.  17).  In  the  second, 
ideal,  case  (Case  4  in  Table  4),  both  ionization  profile  and  magnetic  field  profile  were  uniform  across  a  semi-infinite 
slab  (Fig.  18).  The  load  factor  in  both  cases  was  A=1.15.  (Calculations  with  other  values  of  k  give  qualitatively 
similar  results). 

Results  in  both  cases  show  that  the  MHD  accelerator  increases  the  flow  velocity,  but  due  to  heating,  local  Mach 
numbers  decrease.  Also,  because  of  heating,  mass  capture  and  total  pressure  actually  decrease.  Thus,  the  use  of 
MHD  accelerator  is  only  detrimental.  Table  4  lists  the  key  parameters  in  all  the  cases  with  geometry  designed  for 
M=10. 


c.  MHD  control  at  M<M(iesign*  generator  regime 

As  suggested  recently,  mass  capture  can  be  increased  by  an  on-ramp  MHD  generator  with  a  magnetic  field  that  has 
components  both  parallel  and  orthogonal  to  the  flow.®  Figs.  19-21  illustrate  the  results  of  our  computations  for  three 
versions  of  such  a  device  for  M=6,  h=26  km,  off-design  regime  (the  design  regime  is  M=10,  h=33  km).  The  key 
computed  parameters  are  listed  in  Table  4  as  Cases  5,  6,  and  7.  For  all  three  versions,  the  MHD  region  was  placed 
between  x=22.3  m  and  x=23.3  m,  the  magnetic  field  components  were  B^  =  B^=-2T  (the  negative  sign,  of  course, 
means  that  the  component’s  direction  is  opposite  to  that  of  the  coordinate  axis),  and  the  load  factor  was  ^0.5. 

In  the  first  version  (Fig.  19),  an  electron  beam  power  deposition  of  Qb=10®  W/m®  was  assumed  uniform  through  a 
“slab”  region  that  extends  vertically  down  from  the  ramp  and  does  not  reach  the  cowl  level.  In  the  second  version 
(Fig.  20),  a  Gaussian  distribution  of  beam  power  across  the  “slab”  in  the  vertical  direction  was  assumed.  The  results, 
as  seen  in  Table  4,  show  that  both  versions  give  results  that  are  close  to  each  other  and  quite  negative:  mass  capture, 
total  pressure  coefficient,  and  compression  ratio  all  decrease.  The  negative  performance  is  due  to  Joule  heating  and 
power  introduced  by  the  beam,  resulting  in  entropy  generation  and  thermal  expansion.  The  positive  effect  of  the 

compressing  force  (y  x  =  ~jyB^  is  too  weak  to  overcome  the  thermal  expansion;  no  less  important  is  that  this 

force  does  not  act  on  the  gas  below  the  ionized  “slab”. 
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In  the  third  version  of  MHD  generator  with  flow  compression,  a  uniform  electron  beam  power  deposition  of  Qi,=10* 
W/m^  was  assumed  through  a  “slab”  region  that  extends  vertically  down  to  infinity  (Fig.  21).  In  this  case,  mass 
capture  and  compression  ratio  did  increase  slightly,  although  stagnation  pressure  decreases  in  about  the  same 
proportion  as  the  mass  capture  increase.  Note  that  practical  realization  of  this  version  would  be  very  problematic. 
First,  a  strong  and  reasonably  uniform  magnetic  field  has  to  extend  several  meters  down  form  the  ramp.  Second, 
electron  beams  can  reach  the  intended  MHD  region  only  along  magnetic  field  lines,  so  that  configuration  like  that 
depicted  in  Fig.  22  should  be  created.  In  this  configuration,  ionization  and  heat  addition  would  also  occur  outside  the 
intended  “slab”,  resulting  in  stagnation  pressure  and  mass  capture  losses.  Third,  electric  current  collection  may  be 
problematic. 

As  is  clear  from  comparison  of  the  modestly  successful  third  version  (Fig.  21)  that  has  a  very  long  “slab”  MHD 
region  with  the  versions  of  Figs.  19  and  20,  where  MHD  regions  are  shorter,  most  of  the  ionized  region  in  the  “slab” 
in  Fig.  21  is  useless  or  even  detrimental  for  the  task  of  mass  capture  increase.  It  is  only  the  combined  effect  of 

heating  and  (j  ^sj  =  -J^B^  force  on  the  gas  region  at  the  cowl  lip  level  and  below  it  that  deflects  streamlines  and 

increases  mass  capture. 

IV.  .  Mass  capture  increase  by  off-body  energy  addition:  a  virtual  cowl 

We  propose  a  new  method  of  decreasing  flow  spillage  and  increasing  mass  capture  at  Mach  numbers  lower  than  the 
design  one.  The  essence  of  the  method  is  to  create  a  heated  region  upstream  of  and  somewhat  below  the  cowl  lip 
(Figs.  23  and  24).  The  incoming  flow  would  be  deflected  by  the  elevated-temperature  and/or  elevated  pressure 
region,  causing  an  increased  mass  flow  into  the  inlet.  In  other  words,  a  thermal  funnel,  or  a  virtual  cowl,  would  be 
used  to  increase  mass  capture. 

The  heated  region  may  be  generated  by  supplying  microwave  or  RF  energy  to  a  volume  pre-ionized  by  a  focused 
laser  or  electron  beam.  Other  possible  means  include  plasma  or  hot-air  jets,  and  external  combustion.  Both  off-body 
(Fig.  23)  and  on-lip  (Fig.  24)  variants  can  be  suggested. 

An  advantage  of  this  method  of  inlet  control  as  compared  with  the  MHD  accelerator  method  is  its  relative  simplicity 
and  the  absence  of  problems  with  electron  beams  and  their  transmission  windows,  massive  superconducting 
magnets,  high-voltage  hardware,  and  interelectrode  arcing  due  to  high  Hall  fields  that  characterize  hypersonic  MHD 
devices  in  cold  air.^  '^  Another  advantage  of  the  new  method  is  that  the  air  entering  the  inlet  would  have  experienced 
little  or  no  Joule  heating,  thus  minimizing  stagnation  pressure  losses.  The  virtual  cowl  method  can  have  additional 
advantages  of  increased  lift  and  decreased  drag. 

Preliminary  results  of  two-dimensional  inviscid  modeling  of  mass  capture  increase  by  off-body  heat  addition  are 
shown  on  the  Fig.  25  and  hsted  in  Table  4.  The  freestream  conditions  were  those  at  Mach  6,  h=26  km,  and  the 
geometry  corresponded  to  shock-on-lip  condition  at  Mach  10.  A  heat  source  was  placed  upstream  and  below  the 
inlet.  The  heating  was  distributed  in  a  Gaussian  stretched  ellipsoidal  shape  (Fig.  25).  The  key  parameters  are  listed 
in  Table  4  for  2  cases:  heat  addition  rate  of  Q=5  MW  and  Q=10  MW  per  meter  in  the  transverse  direction  (these 
power  levels  are  equal  to  ~1%  and  ~2%  of  total  enthalpy  entering  the  inlet  in  1  second). 

Mass  capture  in  the  two  cases  increases  by  4-8%.  Interestingly,  both  total  pressure  recovery  and  compression  ratio 
also  improve.  The  reason  for  such  favorable  behavior  of  total  pressure  is  that  the  angle  at  which  gas  enters  the  inlet 
becomes  shallower  due  to  interaction  between  the  deflected  flow  and  the  weak  shock  wave  produced  by  the  heat 
source,  on  one  side,  and  the  ramp  compression  shocks  and  the  flow  downstream  of  the  second  ramp  on  the  other 
side. 

In  the  computed  cases,  the  location,  shape,  and  power  of  the  heat  addition  were  not  optimized.  These  issues  should 
be  investigated  in  future  work. 
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V.  Conclusions 


Among  the  various  scenarios  of  scramjet  inlet  control  investigated  in  this  paper,  two  appear  attractive.  One  is  the  on- 
ramp  MHD  generator  that  would  bring  the  shocks  back  on  lip  at  Mach  numbers  greater  than  the  design  one.  Despite 
the  weight  of  magnets  and  other  components,  overall  weight  savings  may  in  principle  result  from  avoiding  a 
mechanical  system  of  variable  geometry  inlet,  and  from  shorter  forebody  designed  for  lower  Mach  numbers.  The 
MHD  device  is  self-powered,  with  ionization  cost  substantially  lower  than  the  generated  power.  Two  serious  issues 
associated  with  this  method  are:  electron  beam  transmission  windows,  and  losses  of  stagnation  pressure  due  to  Joule 
heating.  Note  that  the  stagnation  pressure  losses  can  be  partially  compensated  with  energy  bypass:  the  generated 
electricity  can  be  use  for  MHD  acceleration  of  the  flow  in  the  combustor  or  downstream  of  it. 

Another  inlet  control  method  that  appears  interesting  is  the  virtual  cowl  concept:  mass  capture  increase  at  Mach 
numbers  lower  than  the  design  one  by  off-body  energy  addition.  While  working  with  fixed-geometry  inlet  and 
avoiding  difficulties  with  electron  beams,  magnets,  and  electrodes  associated  with  MHD  devices,  the  virtual  cowl 
method  not  only  increases  mass  capture,  but  it  also  increases  total  and  static  pressure  at  the  inlet.  The  virtual  cowl 
method  does  require  power  input.  Although  the  required  power  may  be  only  a  small  fraction  of  the  enthalpy  flux  in 
the  inlet,  the  absolute  power,  as  typical  for  hypersonic  regime,  can  be  high.  One  scenario  that  can  be  explored  in 
future  work  is  extraction  of  a  small  fraction  of  flow  enthalpy  with  an  MHD  generator  downstream  of  the  combustor 
and  using  this  power  to  increase  mass  capture  and  compression  at  the  inlet  by  virtual  cowl. 
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Table  1.  Principal  inlet  parameters  in  design  and  off-design  cases  with  and  without  MHD  generator  control.  Design 
conditions:  M=6,  h=26  km,  0.546  atm  freestream  dynamic  pressure.  Off-design  conditions:  M=8,  h=30  km,  0.53 
atm  freestream  dynamic  pressure,  and  M=10,  h=33  km,  0.527  atm  freestream  dynamic  pressure.  '  ,  and  T  are 

the  pressure,  total  pressure  and  translational  temperature,  averaged  over  the  inlet  entrance  cross  section. 


Case 

m, 

kg/(ms) 

H, 

MW/m 

Compression 
ratio:  p! 

Total  pressure 
recovery:  p,lp,^^ 

f  ,K 

Mass 

capture 

M=6,  h=26  km 
design 

721 

49.8 

0.589 

785 

1 

M=8,  h=30  km 
no  MHD 

284.4 

883 

64.9 

0.432 

916.8 

1 

M=8,  h=30  km 
with  MHD 

284.4 

872 

75.6 

0.223 

1235 

1 

M=10,h=33km 
no  MHD 

224.6 

1070 

81.6 

0.305 

1128 

1 

M=10,h=33km 
with  MHD 

224.6 

1030 

108.7 

0.098 

1832 

1 

Table  2.  MHD  generator  characteristics  for  the  cases  of  Table  1.  Electron  beam  current  density  is  y,,  =1  mA/cm^, 
the  beam  relaxation  length  is  =  2.3  m,  and  the  MHD  load  factor  is  ^0.5. 


Case 

,  MW/m 

P-.  g ,  MW/m 

jxB  ’ 

Pj,MW/m 

M=8,  h=30  km 

20.9 

5.42 

41.9 

16.8 

2.4 

M=10,  h=33  km 

42.9 

4.43 

85.8 

29.3 

4.16 

Table  3.  Estimated  displacement  thickness,  99%  velocity  thickness,  and  thermal  boundary  layer  thickness  in  Mach  6 
design  and  Mach  8  off-design  (with  and  without  MHD)  cases. 


Results  at  x=10  m  downstream 
from  vehicle  nose 

Sj. ,  cm 

J„,cm 

M=6;  h=26  km;  design 

5.08 

12.0 

7.01 

M=8;  h=30  km;  no  MHD 

6.08 

13.42 

7.674 

M=8;  h=30  km;  with  MHD 

5.71 

13.46 

10.4 

Table  4.  Performance  of  various  methods  of  inlet  control  at  Mach  number  lower  than  the  design  one.  Design 
conditions:  M=10,  h=33  km.  All  other,  ofF-design,  regimes:  M=6,  h=26  km.  p,p^ ,  and  T  are  the  static  pressure, 
total  pressure,  and  translational  temperature,  averaged  over  the  inlet  cross  section;  Q  is  the  heat  addition  rate  per  unit 
length  in  the  transverse  direction,  with  heat  source  located  upstream  and  below  the  cowl  lip 


Case 

number 

Case  description 

m, 

kg/(ms) 

H, 

MW/m 

Compression 
ratio:  p!  p^ 

Total 

pressure 

recovery; 

7',K 

Mass 

capture 

1 

261 

131.6 

0.31 

1245 

1 

2 

noMHD;Q=0 

274.9 

mi 

46.1 

0.607 

758.6 

0.585 

3 

MHD  accelerator, 
“short  slab”  Gaussian 
ionization  profile 

237.9 

43.0 

0.494 

835.4 

4 

MHD  accelerator, 
“infinite  slab” 
uniform  ionization 

264.9 

1^1111 

47.53 

0.549 

818 

5 

MHD  generator, 
“infinite  slab” 
uniform  ionization 

283.9 

519 

48.81 

0.592 

773.2 

0.593 

6 

MHD  generator, 
“short  slab”  uniform 
ionization  profile 

266.4 

487 

45.1 

0.578 

763 

0.557 

7 

MHD  generator, 
“short  slab”  Gaussian 
ionization  profile 

266.2 

487 

45.2 

0.572 

769 

8 

Virtual  cowl, 
e=5  MW/m 

289.4 

530 

47.9 

0.632 

760.9 

9 

Virtual  cowl, 

0=10  MW/m 

302.5 

554 

50.2 

0.642 

772 

0.632 
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Fig.  1.  Design  forebody  and  inlet  geometry  with  shock-on- 
lip  condition. 


Cowl  lip 

Fig.  2.  Flow  and  shock  geometry  at  Mach  number 
higher  than  the  design  one,  with  no  MHD  control. 


Fig.  3.  Flow  and  shock  geometry  at  Mach  number  lower 
than  the  design  one,  with  no  MHD  control. 


Cowl  lip 


Fig.  4.  Flow  and  shock  geometry  at  Mach  number 
higher  than  the  design  one,  with  MHD  generator 
control. 


Cowl  lip 

Fig.  5.  Flow  and  shock  geometry  at  1  VJIUVII  llUlll  oer  lower 
than  the  design  one,  with  MHD  accelerator  control. 
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Fig.  7.  Comparison  of  electron  beam  power  deposition  profiles  predicted  by  Monte-Carlo  calculations  (MCC)  with 
CYLTRAN,  Gaussian  approximation  with  empirical  relaxation  length  and  Gaussian  approximation 

with  the  given  by  the  MCC. 
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X,  m 


X,  m 


Fig.  9.  Shock  pattern  (upper  plot)  and  flow  streamlines  (lower  plot)  in  M=8,  h=30  km,  off-design  regime,  without 
MHD  control,  for  the  geometry  of  Fig.  8. 
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X,  m 


Fig.  10.  Pressure  contours,  temperature  contours,  flow  streamlines,  and  electron  beam  power  deposition  contours  in 
Mach  8,  h=30  km,  case  with  MHD  generator  control,  for  the  geometry  of  Fig.  8. 
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Fig.  11.  Contour  lines  of  number  densities  of  electrons  and  positive  and  negative  ions,  and  vibrational  temperatures 
contour  lines  in  Mach  8,  h=30  km,  case  with  MHD  generator  control,  for  the  geometry  of  Fig.  8. 
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X,  m 

Fig.  13.  Electron  beam  energy  at  the  injection  point  versus  coordinate  along  the  MHD  generator  channel  in  the  cases 
of  Figs.  10  and  1 1  and  Fig.  12,  with  =2.3  m  =  const. 


z,  cm 

Fig.  14.  Computed  temperature  and  velocity  profiles  in  the  boundary  layer  at  r=10  m  downstream  of  the  vehicle 
nose  in  design  regime  of  Mach  6,  h=26  km  (the  case  of  Fig.  8),  and  off-design  regime  M=8,  h=30  km,  without  (the 
case  of  Fig.  9)  and  with  MHD  (the  case  of  Figs.  10  and  1 1). 
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Fig,  17.  Local  Mach  number  and  horizontal  flow  velocity  contours  in  M=6,  h=26  km,  ofF-design  regime  (geometry 
of  Fig.  15),  with  MHD  accelerator  parameters:  B=4  T,  t=1.15,  andi./f=2.9  m. 
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Fig.  18.  Local  Mach  number  and  horizontal  flow  velocity  contours  in  M=6,  h=26  km,  off-design  regime  (geometry 
of  Fig.  15),  with  MHD  accelerator  at  5=4  T,  ^1.15,  and  with  uniform  electron  beam  power  deposition  of  Qb=10* 
W/m^  in  semi-infinite  “slab”  located  fromA:=8.5  to;t=10.5  m. 
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Fig.  20.  Electron  beam  power  deposition  and  local  Mach  number  contours  in  MHD  generator  with  flow 
compression  case  M=6,  h=26  km,  off-design  regime  (geometry  of  Fig.  15).  Load  factor  A?=0.5;  magnetic  field 
components  Bx=B^=-2  T.  Beam  power  deposition  along  z  is  Gaussian  in  “short  slab”  from  r:=22.3  up  to  r=23.3  m. 


Fig.  21.  Local  Mach  number  contours  in  MHD  generator  with  flow  compression  case  M=6,  h=26  km,  off-design 
regime  (geometry  of  Fig.  15).  Load  factor  k=0.5;  magnetic  field  components  B/^B^—l  T.  Beam  power  deposition  of 
Qb=10®  W/m^  is  uniform  in  semi-infinite  “slab”  fromx=22.3  up  tor=23.3  m. 
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Fig.  24.  Schematic  diagram  of  the  virtual  cowl  concept;  heat  addition  at  the  cowl  lip  increases  mass  capture. 
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Fig.  25.  Mach  number  contour  lines  in  virtual  cowl  case.  Freestream  conditions:  Mach  6,  h=26  km.  Geometry  is  that 
of  Fig.  15.  Heat  addition  rate  of  Q=5  MW  per  1  meter  in  the  transverse  direction  is  distributes  as  a  stretched 
ellipsoidal  Gaussian  function  of  {x,  z). 
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CHAPTER  4 


DEVELOPMENT  OF  A  SUPERSONIC  PLASMA  WIND  TUNNEL 

ABSTRACT 

The  study  of  controlled  plasma  discharges  in  supersonic  air  has  become  an  area  of  increasing  interest. 
These  discharges  may  be  used  for  flow  control,  drag  reduction,  combustion  initiation  and  stabilization,  as  well  as  for 
power  extraction  and  plasma  shielding  applications.  In  order  to  study  the  formation  of  such  discharges  and  their 
impact  on  the  flow,  proper  testing  facilities  are  required.  This  chapter  details  the  work  done  in  the  design, 
construction,  and  testing  of  a  supersonic  plasma  wind  tunnel. 

1.  Introduction 

Over  the  past  decade,  research  into  the  interaction  of  plasma  with  a  gas  flow  has  yielded  some  results 
which  have  not  yet  been  fully  explained’’*.  For  example,  changes  in  shock  wave  strength  and  shock  standoff 
distance  caused  by  the  presence  of  plasma  have  been  observed  experimentally,  but  cannot  be  fully  explained  with 
present  theoretical  models.  A  better  understanding  of  these  interactions,  as  well  as  the  potential  applications  for  the 
controlled  modification  of  shock  structure,  motivates  the  present  work. 

IL  Wind  Tunnel  Design 

The  basic  configuration  of  the  wind  tunnel  is  controlled  by  the  need  to  create  a  plasma  in  the  region  of 
supersonic  flow.  A  schematic  of  the  facility  is  shown  in  Figure  1  below. 

A  plasma  is  created  with  a  50  kW,  1  ms  pulse  of  microwave  radiation  at  2.45  GHz,  which  induces  electric 
breakdown  of  the  air  in  the  region  of  supersonic  flow.  Microwaves  are  introduced  into  the  plenum  through  an  EM 
window  mounted  over  a  port  on  one  wall.  The  direction  of  propagation  is  then  turned  90°  by  an  aluminum  reflector, 
so  the  microwaves  travel  coaxially  with  the  flow  through  the  plenum  chamber  and  nozzle  section.  The  aluminum 
reflector  is  made  with  a  dense  array  of  %”  holes  in  it  to  allow  passage  of  air  while  simultaneously  reflecting  the 
microwaves.  The  breakdown  location  is  controlled  by  two  factors:  a  reduced  static  pressure  in  the  high  speed 
section,  and  an  increased  field  intensity  downstream  of  the  throat. 

Microwave  discharges  occur  more  easily  at  low  pressure  due  to  the  reduced  collision  frequency.  Thus, 
introducing  microwaves  through  the  plenum  where  the  air  is  at  high  pressure  avoids  the  risk  of  parasitic  breakdovm 
near  the  window  housing  and  window  surface. 

As  shovm  in  Figure  1,  the  waveguide  walls  are  tapered  to  create  an  enhancement  of  the  electric  field  as  the 
microwaves  propagate  downstream.  The  nozzle  contours  within  the  taper  are  made  of  dielectric  material  so  as  to 
transmit  the  microwave  radiation.  The  taper  continues  dovra  to  and  slightly  beyond  the  cutoff  dimension  for 
radiation  at  2.45  GHz.  This  generates  a  reflection,  reversing  the  propagation  direction  as  the  microwaves  reach  the 
cutoff  point.  This  taper  is  used  as  one  end  of  a  microwave  resonator.  The  other  end  of  the  resonator  is  formed  by  a 
triple  stub  tuner  on  the  generator  side  of  the  entrance  window.  This  arrangement  allows  for  the  build  up  of  a 
standing  wave  in  the  wind  tunnel.  In  order  to  provide  a  further  field  enhancement  at  the  intended  breakdown 
location,  and  to  provide  for  spatial  localization  of  the  plasma,  two  1  cm  long  metal  pins  are  mounted  opposite  each 
other  on  the  inside  walls  of  the  wind  tunnel.  The  magnetron  is  protected  from  spurious  reflections  by  a  circulator 
and  dummy  load. 

The  flow  of  air  proceeds  from  the  plenum,  through  the  nozzle,  and  into  the  2”x2”  test  section.  The  air  is 
then  discharged  into  a  vacuum  tank.  The  nozzle  is  contoured  to  achieve  Mach  3  flow  in  the  test  section.  Flow 
conditions  in  the  test  section  are  20  Torr  static  pressure,  1 10  K  static  temperature,  and  a  Mach  number  of  3. 

III.  Microwave  Field  Distribution 

The  microwave  field  distribution  was  measured  along  the  centerline  of  the  wind  tunnel  from  the  reflector 
plate  to  the  wind  tunnel  exit  plane.  For  these  measurements,  the  microwave  field  was  generated 
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figure  2:  PLENUM  AND  NOZZLE  CUTAWAY  DRAWING 
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using  a  100  mw  source  operating  in  a  continuous  mode  at  2.45  GHz.  This  source  was  coupled  into  the  wind  tunnel 
with  the  same  triple  stub  tuner  and  window  arrangement  used  for  the  50  kW,  1  ms  pulse  breakdown  experiments. 
The  measurements  were  made  with  a  probe  consisting  of  a  1  cm  antenna  protrading  from  a  length  of  semi-rigid 
coaxial  cable.  This  probe  was  introduced  into  the  wind  tunnel  from  the  downstream  end  and  traversed  through  the 
wind  tunnel  in  increments.  The  2.45  GHz  signal  obtained  with  the  probe  was  attenuated  and  rectified  with  a  crystal 
diode  before  being  recorded  on  a  digital  oscilloscope.  An  average  of  several  such  traverse  measurements  is  shown 
in  Figure  3. 

IV,  Numerical  Simulations 

To  complement  the  measurements  described  above,  as  well  as  to  permit  parametric  studies  of  the  effect  of 
changes  in  the  electrical  properties  of  the  wind  tunnel,  a  commercial  finite  element  code  was  used  to  simulate  the 
electric  field  standing  wave  pattern. 

The  code  used  was  the  High  Frequency  Structure  Simulator  package  from  Hewlett  Packard,  which  was  run 
on  a  UNIX  based  workstation.  The  results  of  several  case  studies,  including  a  comparison  with  the  measurements 
described  above,  are  presented  in  Figures  4  through  8. 

In  the  above  simulations,  incident  wave  phase  is  chosen  to  show  the  standing  wave  pattern  at  its  maximum 
intensity.  For  the  cases  shown  in  Figures  4  and  5  the  region  of  highest  field  intensity  remains  localized  at  the  same 
point  throughout  each  cycle  of  the  incident  wave  propagation.  For  the  cases  shown  in  Figures  6  and  7,  due  to  the 
lossy  dielectric,  deviations  from  a  perfect  standing  wave  pattern  occur.  In  particular,  the  region  of  highest  field 
intensity  changes  dramatically  as  the  incident  wave  goes  through  one  cycle.  The  field  can  be  seen  to  be  strongly 
dependent  on  dielectric  material  properties,  which  has  important  implications  for  the  generation  of  plasma  inside  the 
wind  tuimel. 

V.  Breakdown  Measurements 

Using  the  G-10  nozzle  breakdown  can  occur  in  either  of  two  locations,  depending  upon  conditions.  Under 
static  conditions  below  about  30  Torr,  breakdown  occurs  upstream  of  the  throat,  slightly  below  the  centerline  of  the 
wind  tunnel.  Above  30  Torr  under  static  conditions  breakdown  occurs  downstream  of  the  throat  at  the  intended 
location.  These  tw'o  positions  correspond  with  the  two  peaks  in  field 

intensity  in  the  simulation.  Note  that  imder  flow  conditions,  with  higher  pressure  upstream  of  the  throat  than 
downstream,  the  breakdown  will  also  occur  in  the  downstream  location. 

Breakdown  in  the  wind  tunnel  was  recorded  with  photodiodes  and  high  speed  photographic  equipment. 
The  pressure  was  chosen  for  the  static  case  so  as  to  match  gas  density  with  the  flow  conditions.  Measurements  of 
the  plasma  luminosity  were  made  both  with  flow  and  under  static  conditions.  In  all  cases  the  tests  were  performed 
in  air  with  a  microwave  pulse  duration  of  1  ms.  The  results  indicate  the  initiation  of  the  plasma  within  10  Ds  after 
the  start  of  the  microwave  pulse  without  flow,  and  within  20  Os  after  the  start  of  the  pulse  with  flow.  A  brighter 
plasma  was  observed  with  flow  than  without  flow,  as  shown  in  Figures  9  and  10. 

The  photographs  shown  in  Figure  1 1  were  made  with  a  grated,  intensified  CCD  camera.  The  view  is 
through  the  axial  diagnostic  port  shown  in  Figure  1,  looking  upstream  through  the  wind  tuimel.  Images  were  taken 
of  the  plasma  both  with  the  wind  tuimel  running  and  under  static  conditions  at  60  Torr  and  room  temperature.  The 
pressure  was  chosen  for  static  conditions  to  match  the  density  with  the  Mach  3  flow  for  comparison.  In  all  cases  the 
camera  gate  was  maintained  at  3  Ds.  The  images  show  the  plasma  generated  while  the  wind  tunnel  was  running  to 
be  brighter  and  more  contracted  than  the  plasma  generated  under  static  conditions.  The  location  of  the  plasma 
within  the  tunnel  appears  to  be  stable  for  the  duration  of  the  1  ms  pulse.  Note  that  each  image  represents  a  separate 
instance  of  plasma  generation. 

VT.  Conclusions 

A  Mach  3  wind  tunnel  with  integrated  microwave  plasma  generating  capability  has  been  designed  and 
tested.  The  importance  of  the  electrical  properties  of  the  wind  tunnel  for  microwave  propagation  and  breakdown  has 
been  shown.  Self  induced  plasma  generation  has  been  achieved  with  a  microwave  pulse  length  of  1  ms.  A  stable, 
well  localized  plasma  has  been  maintained  in  the  wind  tunnel  with  flow.  Measurements  of  the  plasma  indicate  a 
brighter  and  more  filamentary  discharge  with  flow  than  without  flow. 


75 


REFERENCES 

1.  W.  Beaulieu,  V.  Bytyrin,  A.  Klimov,  S.  Leonov,  A.  Pashina,  and  B.  Timofeev,  “Plasma  Aerodynamic 
Wind  Tunnel  Tests  with  1/6  Scale  Model  of  Nose  Part  of  F-15,”  3’^“'  Weakly  Ionized  Gases  Workshop, 
Norfolk,  VA,  November  1-5,  1999. 

2.  A.  Klimov,  S.  Leonov,  A.  Pashina,  V.  Skvortsov,  T.  Cain,  and  B.  Timofeev,  “Influence  of  a  Corona 
Discharge  on  the  Supersonic  Drag  of  an  Axisymmetric  Body,”  3'“*  Weakly  Ionized  Gases  Workshop, 
Norfolk,  VA,  November  1-5,  1999. 

3.  S.  Merriman,  I.  Adamovich,  and  J.W.  Rich,  “Studies  of  Oblique  Shock  Waves  in  Weakly  Ionized 

Nonequilibrium  Plasmas,”  3'^‘‘  Weakly  Ionized  Gases  Workshop,  Norfolk,  VA,  November  1-5,  1999. 

4.  R.  Yano,  V.  Contini,  E.  Ploenjes,  P.  Palm,  S.  Merriman,  I.  Adamovich,  W.  Lempert,  V. 
Subramaniam,  and  J.W.  Rich,  “Flow  Visualization  in  a  Supersonic  Nonequilibrium  Plasma  Wind  Tunnel,"  30”' 
Plasmadynamics  and  Lasers  Conference,  Norfolk,  VA,  June,  1999. 

5.  V.  Soloviev,  V.  Krivstov,  A.  Konchakov,  N.D.  Malmuth,  "Mechanisms  of  Shock  Wave  Dispersion 
and  Attenuation  in  Weakly  Ionized  Cold  Discharge  Plasmas,"  Perspectives  of  MHD  and  Plasma 
Technologies  in  Aerospace  Applications  Conference,  IVTAN,  Moscow,  Russia,  March  24-25,1999. 

6.  S.  Macheret,  L.  Martinelli,  and  R.  Miles,  "Shock  Wave  Propagation  in  Non-Uniform  Gases  and 
Plasmas,"  37’*'  AIAA  Aerospace  Sciences  Meeting  and  Exhibit,  Reno,  NV,  January  1 1-14,  1999. 

7.  A.  Yalin,  Y.  lonikh,  R.  Miles,  "Temperature  Measurements  in  Glow  Discharges  with  Ultraviolet 
Filtered  Rayleigh  Scattering,"  30”'  AIAA  Plasmadynamics  and  Lasers  Conference,  Norfolk,  VA,  June  28-July 
1,  1999. 

8.  P.  Harsha  and  E.  Gurjanov,  "AJAX:  New  Directions  in  Hypersonic  Technology,"  8”'  Aerospace  Planes 
and  Hypersonic  Technologies  Meeting,  Norfolk,  VA,  April,  1996. 

9.  A.D.  MacDonald,  Microwave  Breakdown  in  Gases,  Wiley,  1966. 

10.  R.  Collin,  Foundations  for  Microwave  Engineering,  McGraw-Hill,  1966. 

11.  G.  Matthaei,  L.  Young,  E.  Jones,  Microwave  Filters,  Impedance  Matching  Networks,  and 
Coupling  Structures,  McGraw-Hill,  1964. 


76 


normalized  signal 


Relative  Field  Strength 


FIGURE  3:  MEASURED  ELECTRIC  FIELD  ALONG  WIND  TUNNEL  CENTERLINE 


77 


t|E| 


FIGURE  4:  ELECTRIC  FIELD  MAGNITUDE  WITH  AIR  FILLED  WAVEGUIDE  (NO  DIELECTRIC) 


FIGURE  5;  ELECTRIC  FIELD  MAGNITUDE  WITH  TEFLON  NOZZLE 
Relative  permittivity  =  2.08  Dielectric  loss  tangent  =  0.00037 
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FIGURE  6:  ELECTRIC  FIELD  MAGNITUDE  WITH  G-10  NOZZLES 
Relative  permittivity  =  5.0  Dielectric  loss  tangent  =  0.019 
Incident  wave  phase  =  157“ 


LE  7:  ELECTRIC  FIELD  MAGNITUDE  WITH  G-10  NOZZLl 
Relative  permittivity  =  5.0  Dielectric  loss  tangent  =  0.019 
Incident  wave  phase  =  70“ 
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FIGURE  8:  COMPARISON  BETWEEN  SIMULATION  AND  MEASUREMENT  USING  G-IO 
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FIGURE  9:  PLASMA  LUMINOSITY  WITHOUT  FLOW  (60  TORE,  300  K) 


Microwave  Power 


Luminosity 


FIGURE  10:  PLASMA  LUMINOSITY  WITH  FL 


>  800  Jis 

FIGURE  1 1 :  COMPARISON  OF  PLASMA  WITH  AND  WITHOUT  FLOW 
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CHAPTER  5 


STUDIES  OF  ELECTRON  BEAM-INITIATED,  MICROWAVE  DRIVEN  PLASMA  DISCHARGES  FOR 

AERODYNAMIC  CONTROL 

ABSTRACT 

Control  of  supersonic  vehicles  by  means  of  heat  addition  is  possible  using  off  body  heat  addition  by 
microwave  breakdown  in  air.  The  resulting  changes  in  temperature  and  shock  stracture  alter  the  pressure 
distribution  around  a  body  downstream  of  the  heat  addition  region.  In  the  present  work,  a  controlled  microwave 
driven  discharge  in  a  small  scale  Mach  3  wind  tunnel  is  studied  to  investigate  its  potential  use  as  a  means  of  vehicle 
control.  The  location  of  the  discharge  is  controlled  by  the  microwave  power,  the  geometry  of  the  wind  tunnel,  and  a 
low  power  electron  beam  injected  into  the  wind  tunnel.  TTie  electron  beam  was  found  to  reduce  the  microwave 
power  required  for  breakdown,  and  to  provide  a  means  of  modulating  the  plasma  in  time. 

L  Introduction 

Localized  heat  addition  upstream  of  a  moving  body  can  affect  the  lift,  drag,  and  moments  about  the  body 
by  altering  the  pressure  distribution  over  the  surface.  This  has  several  potential  advantages  over  conventional 
aerodynamic  controls.  Unlike  conventional  control  surfaces,  this  type  of  control  may  be  implemented  without  the 
use  of  moving  parts,  and  thus  has  reduced  time  response  and  increased  frequency  range  as  compared  to  conventional 
control  surfaces.  This  type  of  control  is  also  inherently  more  flexible.  The  pressure  distribution  may  be  altered 
anywhere  on  the  body,  rather  than  along  control  surfaces  at  fixed  points  on  the  body. 

The  heat  addition  profile  creates  a  corresponding  perturbation  in  temperature,  Mach  number,  and  stagnation 
pressure  encountered  by  the  body.  Shock  waves  produced  by  the  expanding  gas  in  the  heat  addition  region  change 
the  flow  properties  in  adjacent  regions,  and  can  also  be  made  to  impinge  on  the  body  itself  The  exact  manner  in 
which  the  flow  is  perturbed  is  determined  by  the  size,  heat  deposition  profile,  and  spatial  location  of  the  heat 
addition  region  relative  to  the  body,  and  also  by  the  way  in  which  these  parameters  are  varied  in  time.  CFD 
simulations  showing  detailed  flow  structure  and  optimization  of  the  heat  deposition  region  for  flow  around  a  conical 
model,  similar  to  the  model  used  in  the  present  experiments,  can  be  found  in  reference  2. 

The  present  work  involved  experiments  performed  in  a  Mach  3  wind  tunnel  with  heat  deposited  in  the  test 
section  by  a  microwave  driven  discharge.  The  control  of  this  discharge  is  accomplished  by  adjustment  of  the 
microwave  power,  the  wind  tunnel  geometry,  and  modulation  of  a  20  kV  electron  beam  injected  into  the  wind 
tunnel. 

H.  Apparatus 

The  wind  tunnel  used  in  this  work  is  of  the  indraft  type,  drawing  air  at  atmospheric  conditions  into  the 
plenum  and  exhausting  to  a  vacuum  line  maintained  at  150  torr.  This  facility  has  been  described  in  some  detail  in  a 
previous  publication'.  The  tunnel  is  comprised  of  a  plenum  section,  a  rectangular  nozzle  with  a  Mach  3  expansion 
ratio,  a  2.15”x2.166”  rectangular  test  section,  and  a  5.5°  half  angle  diffuser.  An  air  ejector  pump  connected  through 
a  vacuum  line  to  the  diffuser  exit  maintains  the  overall  pressure  drop  through  the  wind  tunnel.  A  schematic  of  the 
wind  tunnel  is  shown  in  Figure  1 . 

Microwaves  are  used  to  provide  local  energy  addition  in  the  tunnel  by  creating  a  microwave  driven 
discharge  in  the  test  section.  The  microwaves  are  generated  externally  by  a  6  kW,  2.45  GHz  magnetron  and  are 
introduced  via  a  waveguide  into  the  plenum  of  the  wind  tunnel  through  a  port  in  the  top  wall.  The  waveguide  has  a 
4.3”x2.15”  cross  section  and  transmits  the  microwaves  in  the  TEjo  mode.  The  plenum  is  made  of  aluminum  and  is 
of  the  same  dimensions  as  the  waveguide,  allowing  the  microwaves  to  propagate  in  the  same  TEio  mode.  A  three 
stub  tuner  near  the  waveguide  entrance  is  used  to  match  the  impedance  of  the  waveguide  to  that  of  the  tunnel, 
minimizing  reflections  back  toward  the  source.  Upstream  of  the  microwave  port,  a  reflector  plate  perforated  with  an 
array  of  !4”  holes  reflects  microwaves  in  the  downstream  direction  and  prevents  microwaves  from  leaking  out  of  the 
wind  tunnel  upstream  of  the  plenum.  The  holes  in  the  reflector  plate  serve  as  an  air  inlet  to  the  plenum,  but  do  not 
allow  the  passage  of  microwaves  since  the  diameter  of  the  holes  is  much  smaller  than  the  cutoff  dimension  for  2.45 
GHz  radiation.  Adding  microwaves  through  the  plenum  rather  than  directly  to  the  test  section  aids  in  avoiding 
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parasitic  breakdown  around  the  entrance  port  because  the  higher  pressure  in  the  plenum  relative  to  that  in  the  test 
section  suppresses  breakdown  in  the  plenum. 


Perforated  reflector 
plate 


Figure  1:  Wind  tunnel  schematic 

Two  identical  contoured  blocks  of  virgin  electrical  grade  Teflon,  along  with  two  flat  aluminum  walls,  form 
the  nozzle.  Teflon  was  chosen  as  the  nozzle  material  due  to  its  excellent  electrical  characteristics,  being  nearly 
transparent  to  microwaves.  The  Teflon  is  mounted  on  angled  aluminum  plates,  which  form  a  taper  toward  the  test 
section.  This  taper  is  extended  to  be  slightly  beyond  the  microwave  cutoff  dimension.  As  the  microwaves  approach 
cutoff,  the  group  velocity  approaches  zero  and  the  electric  field  intensity  increases.  The  cutoff  taper  also  reflects  the 
incident  microwave  power  back  upstream.  With  a  properly  adjusted  three  stub  timer,  the  wind  tuimel  can  be  made 
to  act  as  a  microwave  cavity  resonator.  A  standing  wave  is  formed  between  the  tuner  and  the  cutoff  taper.  Incident 
microwave  energy  is  stored  in  the  standing  wave,  and  the  electric  field  strength  inside  the  tunnel  becomes  much 
stronger  than  the  field  in  the  waveguide,  with  the  field  near  the  cutoff  point  being  the  strongest.  This  is  the  region 
where  the  discharge  forms. 

The  test  section  is  constracted  of  aluminum  and  has  a  2.15”x2.166”  rectangular  cross  section,  with  1.5” 
diameter  quartz  windows  in  the  side  walls.  Aluminum  tubes  smaller  than  the  cutoff  dimension  are  mounted  over  the 
windows  to  prevent  microwave  leakage  but  still  allow  optical  access.  The  cutoff  taper,  containing  the  region  of 
highest  electric  field  strength,  is  just  upstream  of  the  windows  in  the  test  section. 

The  diffuser  used  downstream  of  the  test  section  is  made  up  of  two  parts.  The  first  part  is  a  6”  long 
aluminum  shock  duct  of  the  same  cross  section  as  the  test  section.  The  second  part  transitions  along  the  axial 
direction  gradually  from  the  rectangular  test  section  cross  section  to  the  circular  vacuum  line  cross  section  and 
expands  through  a  5.5°  half  angle  to  form  a  smooth  transition  to  the  6”  diameter  vacuum  line.  It  is  made  of  stainless 
steel  sheet  metal,  and  contains  a  1”  diameter  stainless  steel  sting  used  for  mounting  models  in  the  wind  tunnel.  The 
sting,  which  also  extends  through  the  shock  duct,  reduces  the  area  for  the  flow.  Together  with  the  attached  model, 
this  arrangement  produces  a  converging  area  for  the  flow  at  the  entrance  to  the  shock  duct.  The  shock  waves  from 
the  model  become  the  initial  stage  of  compression  for  the  diffusfer,  improving  flie  pressure  recovery.  The  sting  also 
acts  as  a  central  conductor,  so  microwaves  that  leak  beyond  the  cutoff  portion  of  the  test  section  can  propagate  along 
the  sting  into  the  diffuser.  A  close  fitting  section  of  aluminum  honeycomb  near  the  downstream  end  of  the  diffuser 
prevents  microwaves  from  passing  into  the  vacuum  line,  while  still  allowing  air  to  flow  through. 

The  electron  beam  assembly  is  located  adjacent  to  the  wind  tuimel.  The  beam  is  generated  by  a  small 
electron  gun,  similar  to  the  gun  in  a  television  picture  tube,  housed  inside  a  vacuum  line.  A  custom  power  supply 
allows  for  the  independent  setting  of  all  grid  and  filament  voltages  on  the  gun,  and  for  switching  the  beam  on  and 
off  at  various  frequencies.  The  electron  beam  is  directed  out  of  the  vacuum  line  via  two  pairs  of  steering  magnets 
through  a  thin  silicon  nitride  window,  and  then  into  the  wind  tunnel.  The  beam  is  directed  parallel  to  the  electric 
field  of  microwaves  and  perpendicular  to  the  flow  velocity.  However,  the  beam  diverges  through  a  large  angle  after 
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passing  through  the  window  due  to  scattering  collisions.  A  simplified  schematic  of  the  electron  beam  assembly  is 
shown  below  in  Figure  2. 

The  electron  beam  window  is  made  from  a  500nm  thick  film  of  silicon  nitride  coated  on  a  silicon  wafer. 
An  etching  process  is  performed  to  remove  a  1x1  mm  square  of  silicon  from  the  wafer,  leaving  only  the  thin  silicon 
nitride  film.  This  film  can  hold  one  atmosphere  of  pressure  across  its  surface  at  room  temperature,  although  heating 
caused  by  the  electron  beam  reduces  its  strength.  Thus  the  required  pressure  that  must  be  held  across  the  window 
limits  the  amount  of  electron  beam  current  which  can  be  transmitted  through  it.  The  window  is  recessed  from  the 
inner  wall  of  the  wind  tunnel  by  one  tenth  of  an  inch  to  provide  some  protection  from  aerodynamic  forces. 


Flow 


Figure  2:  Electron  beam  assembly 


III.  Experiment 

The  experiments  performed  in  the  wind  tunnel  were  all  done  using  one  of  two  models  in  the  test  section:  a 
1"  diameter  15°  half  angle  steel  cone  and  a  similar  cone  made  of  Teflon.  It  was  found  that  with  the  steel  cone,  a 
visible  corona  had  a  tendency  to  form  at  the  tip,  as  can  be  seen  in  Figure  3.  The  steel  cone  was  replaced  with  a 
Teflon  cone  of  the  same  dimensions  to  prevent  the  corona  from  forming  and  to  determine  how  that  might  change  the 
electron  beam  control  of  the  plasma.  The  majority  of  experiments  involved  microwave  power  levels  near  the 
breakdown  threshold,  ranging  from  approximately  900  to  2400  watts  of  transmitted  power.  Once  initiated,  the 
plasma  absorbed  about  60%  of  the  transmitted  power.  Data  was  collected  using  an  RF  crystal  mounted  on  the 
reflected  power  port  of  the  isolator  to  measure  reflected  power,  a  photodiode  to  determine  the  light  emission 
intensity,  and  a  fast  shutter  CCD  camera  to  obtain  images  of  the  plasma. 

Three  distinct  types  of  discharge  control  experiments  were  undertaken:  the  first  type  was  simply  to 
establish  the  breakdown  threshold  both  with  and  without  the  electron  beam,  the  second  was  to  investigate  the 
response  of  the  plasma  to  frequency  modulation  of  the  electron  beam,  and  the  third  was  an  attempt  at  using  the 
electron  beam  to  move  the  plasma  to  a  new  location. 

For  the  breakdown  tests,  the  electron  beam  was  first  turned  off  and  microwave  power  reduced  from  a  value 
well  above  threshold  to  a  value  somewhat  below  threshold  for  breakdown.  These  tests  were  then  repeated  with  the 
electron  beam  set  to  various  current  levels.  The  electron  beam  was  maintained  at  a  constant  20  kV  for  these  tests. 
Due  to  a  power  ripple  present  in  the  output  of  the  magnetron,  fluctuations  in  the  plasma  at  60  Hz  and  higher 
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frequencies  were  present  in  most  of  these  experiments.  This  is  particularly  apparent  near  the  breakdown  threshold, 
where  the  plasma  is  alternately  initiated  and  extinguished  as  the  magnetron  power  fluctuates.  This  is  also  where  the 
electron  beam  had  its  greatest  effect  on  the  plasma,  so  the  data  presented  in  the  results  reflect  this  unintentional 
fluctuation. 

Frequency  modulation  of  the  electron  beam  was  performed  at  microwave  power  levels  for  which  the  beam 
was  required  to  initiate  breakdown,  as  determined  by  the  breakdovm  tests  described  above.  Various  frequencies 
were  used,  and  the  response  of  the  plasma  was  determined  by  reflected  microwave  power  and  light  emission 
measurements. 

The  final  set  of  experiments  to  control  the  discharge  location  were  done  by  inserting  a  metal  pin  into  the 
wind  tunnel  on  the  opposite  side  from  the  electron  beam.  The  pin  acts  to  initiate  breakdown  on  its  side  of  the  wind 
tunnel.  Due  of  its  strong  effect  on  initiation  as  compared  to  the  electron  beam,  the  discharge  showed  a  high  degree 
of  sensitivity  to  the  pin  insertion  depth.  Because  of  this,  the  pin  was  inserted  a  variable  depth  into  the  wind  tunnel 
by  means  of  a  translation  stage  which  could  be  adjusted  in  small  increments.  Experiments  were  first  conducted  with 
the  electron  beam  off,  and  the  pin  inserted  incrementally  until  the  plasma  began  to  form  on  the  side  with  the  pin 
only.  This  was  then  repeated  with  the  electron  beam  on  to  move  the  plasma  over  to  the  other  side  of  the  tunnel. 


Figure  3:  Plasma  over  a  steel  cone  in  Mach  3  flow 
IV.  Results 

The  breakdown  threshold  was  significantly  reduced  by  the  presence  of  the  electron  beam.  Self  initiated 
breakdown  occurred  at  about  1700  watts  of  transmitted  microwave  power.  Using  the  electron  beam,  the  power  level 
required  for  breakdown  dropped  from  1700  watts  to  1100  watts,  a  reduction  of  35%.  This  change  in  breakdown 
threshold  is  not  a  strong 

function  of  the  beam  current,  at  least  in  the  range  tested;  there  was  no  significant  change  in  the  threshold  at  beam 
currents  ranging  from  0.25  DA  to  40  OA.  The  lower  limit  of  this  range  represents  the  practical  lower  limit  for 
which  beam  current  can  be  controlled  with  the  present  facility,  and  the  upper  limit  is  the  maximum  current 
obtainable  without  damage  to  the  electron  gun.  No  changes  in  the  plasma  due  to  the  electron  beam  were  observed  at 
microwave  pow'er  levels  above  that  required  for  self  initiated  breakdown. 

At  power  levels  between  1100  and  1700  watts,  the  electron  beam  was  used  to  modulate  the  plasma  in  time. 
The  electron  beam  was  cycled  on  and  off  at  known  frequencies,  and  the  plasma  luminosity  was  recorded  with  a 
photodiode.  The  plots  presented  in  Figures  4a,  4b,  and  4c  show  the  results  of  these  tests  at  5,  10,  and  20  Hz.  At  20 
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Hz,  the  ripple  in  the  microwave  power  noticeably  interferes  with  the  modulation,  as  an  additional  60  Hz  frequency 
from  the  microwave  power  supply  is  overlaid  on  top  of  the  electron  beam  modulation  frequency. 


Plasma  Luminosity 
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Figure  4a:  Electron  beam  modulated  at  5  Hz 
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Figure  4c:  Electron  beam  modulated  at  20  Hz 
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Figure  5:  Transmitted  microwave  power 

The  power  supply  ripple  is  characterized  in  Figure  5,  showing  a  relatively  large  power  fluctuation  of  ±20%. 
Four  frequencies  at  60  Hz,  360  Hz,  2.3  kHz,  and  27.6  kHz  were  identified  from  the  data.  A  portion  of  a  typical  plot 
of  plasma  luminosity  is  shown  in  Figure  6,  with  fluctuations  at  60  Hz  and  360  Hz  also  apparent  in  the  plasma. 
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Figure  6:  Plasma  fluctuations  due  to  power  supply  ripple 
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The  results  of  third  type  of  experiment  to  control  discharge  location  with  the  electron  beam  were 
inconclusive.  As  described  in  the  previous  section,  a  pin  was  placed  in  the  wind  tunnel  along  the  opposite 
wall  from  the  electron  beam  window.  This  provided  a  means  to  initiate  the  discharge  on  the  pin  side  of  the 
wind  tunnel.  The  electron  beam  was  then  turned  on  in  an  attempt  to  extend  the  discharge  over  to  the  other 
wall.  However,  whenever  the  pin  was  inserted  far  enough  to  reliably  initiate  the  discharge,  the  electron 
beam  had  no  effect.  Conversely,  when  the  pin  was  retracted  from  this  point,  the  discharge  had  a  tendency 
to  jump  between  the  two  side  walls  whether  the  electron  beam  was  present  or  not.  This  could  be  due  to 
vibration  in  the  pin  while  the  wind  tunnel  is  running. 

All  experiments  were  repeated  for  both  a  steel  cone  in  the  test  section  and  a  Teflon  cone  of  the 
same  dimensions.  The  results  were  found  to  be  essentially  the  same  for  both  materials,  despite  the 
formation  of  a  corona  on  the  tip  of  the  steel  cone,  which  did  not  occur  with  the  Teflon  cone. 

V.  Conclusons 

A  stable  microwave  driven  discharge  was  maintained  in  a  small  scale  wind  tuimel  against  a  Mach 
3  flow.  A  20  kV  electron  beam  injected  into  the  wind  tunnel  reduced  the  breakdown  threshold  by  35% 
from  1700  to  1 100  watts.  Within  this  range  of  microwave  power  the  discharge  was  controlled  in  time  by 
modulating  the  electron  beam.  The  discharge  behavior  was  not  a  strong  function  of  beam  current  between 
0.25  DA  and  40  DA,  and  the  electron  beam  did  not  alter  the  discharge  behavior  at  microwave  power  levels 
above  that  required  for  self  initiated  breakdown. 

REFERENCES 


1 .  B.  McAndrew,  R.  Murray,  M.  Shneider,  J.  Fox,  J.  Kline,  P.  Efthimion,  and  R.  Miles,  “Comparison 
of  Numerical  and  Experimental  Results  from  Localized,  Microwave  Driven  Plasma  Energy  Addition  into  a 
Mach  3  Flow”,  AIAA  paper  #2001-3061 . 

2.  I.  Girgis,  M.  Shneider,  S.  Macheret,  G.  Brown,  and  R.  Miles,  “Creation  of  Steering  Moments  in 
Supersonic  Flow  by  Off-Axis  Plasma  Heat  Addition”,  AIAA  paper  #2002-0129. 

3.  Riddle,  L.,  “Electron  Beams  for  Control  of  Microwave  Generated  Plasma”,  Samoff  Corporation 
technical  report,  March  3,  2000. 

4.  Ghose,  R.  N.,  Microwave  Circuit  Theory  and  Analysis.  McGraw-Hill  Book  Co.,  Inc.,  1963. 

5.  B.  McAndrew,  P.  Barker,  and  R.  Miles,  “Development  of  a  Supersonic  Plasma  Wind  Tunnel”, 
AIAA  paper  #2000-0533. 


CHAPTER  6 

EFFICIENT  PLASMA  GENERATION  BY  REPETITIVE  HIGH-VOLTAGE  NANOSECOND  PULSES 

FOR  MHD  FLOW  CONTROL 


I.  INTROUDUCTION 

An  important  factor  in  MHD  flow  control  is  that  in  the  flight  regime  of  interest,  static  temperature  is  too  low  for 
thermal  ionization  of  air.  Therefore,  nonequilibrium  ionization  in  MHD  devices  must  be  done,  and,  because  of  high 
recombination  rates,  the  ionization  must  be  sustained  throughout  the  MHD  region.  The  energy  cost  of  ionization 
imposes  very  rigid  constraints  on  the  choice  of  ionization  methods:  only  the  most  energetically  efficient  ionization 
method  can  be  used  in  these  devices.  Our  earlier  work  [3-10]  showed  that  electron  beams  injected  into  the  gas  along 
magnetic  field  lines  constitute  the  most  efficient  way  of  creating  nonequilibrium  ionization. 

An  alternative  method  is  to  generate  high-energy  ionizing  electrons  inside  the  plasma  by  applying  a  very  high 
(super-breakdown)  voltage  for  a  very  short  (nanosecond-scale)  time,  and  to  repeat  such  pulses  at  a  rate  matching  the 
recombination  rate  [10-12].  Theoretical  analysis  [10-12]  demonstrated  that  the  repetitive-pulse  method  is  much 
more  efficient  than  is  ionization  by  a  steady-state  DC  or  RF  field,  but  less  efficient  than  the  ionization  be  electron 
beams.  The  repetitive-pulse  method  is,  however,  free  from  window  problems  associated  with  electron  beams.  Fig.  1 
shows  that  energy  cost  per  newly  generated  electron  as  a  function  of  the  reduced  electric  field,  and  Fig.  2  shows 
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theoretical  prediction  of  power  budget  required  to  sustain  an  average  electron  density  of  lO’^  cm'^  in  room- 
temperature  air  at  a  pressure  of  10  Torr,  and  of  the  required  pulse  repetition  rate  [10-12], 

A  custom-developed  pulser  that  produces  1-3  nanosecond  pulses  with  30  kV  peak  voltage  at  up  to  100  kHz 
repetition  rate,  with  average  power  of  500-800  W,  has  been  purchased  and  installed.  Fig.  3  is  a  photograph  of  the 
plasma  produced  with  this  pulser  between  two  flat-plate  electrodes  in  air  at  1 0  Torr  at  50  kHz  repetition  rate.  This 
plasma  generator  has  been  integrated  with  a  Mach  4  nozzle  and  2.5  inch  working  section  of  a  6.5  Tesla 
superconducting  magnet  purchased  earlier.  Experiments  on  MHD  effects  in  supersonic  cold  air  with  nonequilibrium 
ionization  have  begun. 


DC  or  RF  plasmas 


Fig.  1.  Energy  cost,  in  eV  per  newly  produced  electron  in  weakly  ionized  air  plasmas  versus  ratio  of  electric  field 
strength  E  to  the  gas  number  density  N.  Values  of  E/N  corresponding  to  bulk  DC  and  RF  plasmas  and  to  the  cathode 
sheath  of  glow  discharges  are  indicated  by  vertical  arrows.  The  energy  cost  of  ionization  by  high-energy  electron 
beams  is  indicated  by  a  horizontal  arrow. 
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E/N,  Vcm 


E/N,  V*cm‘ 


Fig.  2.  Power  input  (upper  plot)  and  pulse  repetition  rate  (lower  plot)  required  to  sustain  a  time-averaged  electron 
number  density  of  lO’^  cm'^  in  air  at  10  Torr,  300  K,  by  repetitive  pulses,  as  functions  of  E/N  at  various  pulse 

lengths  r, .  Power  inputs  required  to  sustain  the  same  electron  density  at  steady  state  by  DC  discharge’  and  high- 
energy  electron  beams  are  indicated  in  the  upper  plot. 


Fig.  8.  Plasma  produced  between  flat-plate  electrodes  in  air  at  10  Torr  by  1  ns,  30  kV  pulses  at  50  kHz  repetitic 
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CHAPTER  7 


SEMICLASSICAL  MODELING  OF  STATE-SPECIFIC  DISSOCIATION 
RATES  IN  DIATOMIC  GASES 


I.  Introduction 

The  vibration-dissociation  coupling  is  of  key  importance  to  hypersonic  reentry,  significantly  affecting 
aerodynamics,  radiative  and  convective  heat  fluxes,  and  spectral  signatures  of  vehicles  flying  at  sub-orbital  to  super¬ 
orbital  velocities  in  rarefied  atmosphere.  Over  the  last  four  decades,  high-temperature  dissociation  and  vibrational 
relaxation  was  a  subject  of  a  number  of  experimental  and  theoretical  studies,  reviewed  by  Park,'  and,  more  recently, 
by  Macheret.^  A  number  of  models,  such  as  those  suggested  by  Marrone  and  Treanor,^  Losev,^  and  Park,'’  ®  have 
been  used  with  some  success  to  describe  experimental  data  on  nonequilibrium  dissociation.  Those  models,  however, 
are  empirical  or  semiempirical,  using  adjustable  parameters  found  from  a  fit  to  a  limited  number  of  experimental 
data,  which  caimot  guarantee  their  validity  outside  of  the  calibration  range. 

Most  models  of  vibrational  relaxation,  until  recently,  relied  upon  a  first-order  perturbation  (FOP)  method,^  such 
as  that  used  in  SSH  theory,^’  ®  with  some  modifications. ^  Despite  a  reasonable  agreement  of  the  calculated  overall 
relaxation  times  with  experiments,  FOP  methods  caimot  pretend  to  adequately  describe  state-specific  relaxation 
rates  at  high  collision  energies  and  temperatures  typical  of  hypersonic  bow  shocks. 

Recently,  the  authors  and  their  colleagues  have  developed  new  theories  of  high-temperature  relaxation'®'''*  and 
dissociation^’  '®,  especially  well  suited  for  high-temperature  environment.  The  principal  idea  of  these  theories  is 
that  at  high  energy  of  colliding  particles  a  molecule  can  jump  over  many  vibrational  quantum  states.  These  jumps 
can  be  both  bound-bound  (vibrational  relaxation)  and  bound-free  (dissociation),  and  they  must  be  correctly 
accounted  for  even  in  the  first  approximation.  For  the  relaxation  problem,  correct  and  convenient  framework  is 
provided  by  the  Forced  Harmonic  Oscillator  (FHO)  model,'’"^®  which  naturally  and  consistently  includes  multi¬ 
quantum  jumps.  With  relatively  simple  corrections  for  anharmonicity  and  enforcement  of  detailed  balance,  this 
model  has  been  demonstrated  to  agree  very  well  with  both  experimental  data  and  sophisticated  state-of-the-art 
quantum  computations.  The  FHO  approach  was  also  used  to  model  dynamics  of  populations  of  vibrational  levels 
during  relaxation,  and  to  resolve  a  long-debated  “bottleneck”  problem.'^  Most  recently,  the  approach  was  extended 
to  explicitly  account  for  three-dimensional  motion  of  colliding  particles  and  molecular  rotation,  resulting  in  a  fully 
analytical,  without  adjustable  parameters,  non-perturbative  theory  of  multi-quantum  vibrational  energy  transfer.'^’  '“* 
The  important  role  of  rotation  and  non-collinear  collisions  was  clearly  demonstrated,  and  the  calculated  rates  agree 
well  with  experiments  as  well  as  with  close-coupled  quantum  and  semiclassical  computations  using  full  potential 
energy  surfaces  for  Na-He,^'"^^  N2-N2,^‘'  and  02-02.^^'^^ 

The  recently  developed  dissociation  model^'  also  argued  that  because  of  availability  of  high  kinetic  energies 
of  particles  in  the  post-shock  environment  molecules  can  dissociate  even  from  low  vibrational  states  in  a  single 
collision.  The  model  used  classical  impulsive  (sudden)  approximation  to  calculate  both  dissociation  thresholds  and 
probabilities  to  find  the  colliding  pair  near  its  optimum  configuration.  Vibrational  state-specific  and  rotational 
energy  dependent  dissociation  rates  have  been  derived  in  closed  analytical  form.  However  successful,  the  model^’ 
is  not  free  from  problems.  The  impulsive  approximation  should  certainly  be  correct  as  a  high-energy  asymptotic, 
but  it  is  not  clear  if  a  real  dissociating  system  is  close  to  this  limit.  Additionally,  the  calculated  probability  factors'* 
turned  out  to  have  singularities  at  certain  energies,  which  calls  for  more  sophisticated  methods  of  evaluating  the 
probabilities. 

In  the  present  chapter,  we  suggest  a  theory  of  state-specific  dissociation  rates  of  simple  diatomic  molecules.  The 
approach  is  essentially  a  modified  FHO  model  that  allows  to  calculate  probabilities  of  jumps  across  many  quantum 
states,  if  an  energy  transferred  to  a  classical  initially  non-vibrating  oscillator  can  be  determined.  The  latter  value  is 
derived  for  three-dimensional  collisions  of  rotating  molecules  in  two  cases:  free-rotation  approximation  and 
impulsive  (sudden)  limit.  As  a  result,  both  state-specific  and  thermally-averaged  dissociation  rates  are  calculated, 
and  roles  of  various  degrees  of  freedom  can  be  quantified. 

H.  The  FHO-based  dissociation  model 

Consider  collisions  of  diatomic  molecules  AB  with  atoms  M  in  the  course  of  which  the  molecule  changes  its 
vibrational  quantum  state  from  i  to  /: 
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AB{i)+M^AB{f)  +  M 


(1) 


If  the  molecule  is  modeled  as  a  harmonic  oscillator,  and  the  collision’s  effect  is  treated  as  exerting  a  perturbing  force 
on  the  oscillator,  then  the  probability  of  the  process  (1)  is:’’ 


(-ly 


[[i-r)\{f-r)\Q' 


(2) 


where  n  =  min(/,  /) ,  and  Q  =  AE  /  ha  is  the  dimensionless  energy  (i.e.  the  average  number  of  quanta)  transferred 

to  the  initially  non-vibrating  classical  oscillator  in  an  AB-M  collision. 

The  Forced  Harmonic  Oscillator  (FHO)  method  is  a  consistent  and  convenient  scaling  tool,  allowing  close- 
coupled  calculations  of  transitions  across  an  arbitrary  number  of  vibrational  levels.  This  method  reduces  the 
complexity  of  the  quantum  problem  to  classical  calculations  of  collisional  energy  transfer  to  an  initially  non¬ 
vibrating  oscillator.  The  principal  drawbacks  of  the  FHO  method  are;  i)  anharmonicity  of  real  molecules-oscillators, 
and  ii)  neglect  of  the  reverse  effect  of  the  oscillator  on  the  perturbing  classical  trajectory.  The  simplest  way  to 
account  for  the  anharmonicity  is  to  correct  only  frequency,  and  not  the  wavefunctions:  for  each  transition  i  — >  /  , 

the  quantity  « =  |£; ,  where  andE^  are  the  energies  of  the  respective  levels,  is  used  as  the 

effective  oscillator  frequency.  To  account  for  the  oscillator  feedback  effect  on  the  classical  trajectory,  a 
s3mmetrization  procedure  is  commonly  used  whereby  the  effective  initial  velocity  of  the  colliding  pair  is  taken  as  an 
arithmetic  average  between  the  actual  velocity  and  that  calculated  with  the  kinetic  energy  reduced  by  the  oscillator 
energy  increase,  E.  -  E^  These  two  simple  corrections  have  proved  to  be  successful  in  bringing  the  FHO- 

calculated  vibrational  energy  transfer  rates'*’ "’  into  an  excellent  agreement  with  more  sophisticated  quantum 
calculations^'’  ’’’ 

To  model  dissociation,  we  need  to  calculate  probabilities  of  bound-free  transitions.  The  original  FHO 
theory,  strictly  speaking,  cannot  describe  such  transitions.  First  of  all,  the  discrete  energy  spectrum  of  harmonic 
oscillator  extends  to  infinity.  Secondly,  even  for  a  tmncated  oscillator,  calculating  probabilities  of  bound-free 
transitions  is  difficult  and  is  not  done  directly  by  the  FHO  theory.  In  this  paper,  we  will  make  use  of  the  fact  that, 
when  molecular  vibrations  are  modeled  by  a  Morse  oscillator  with  the  lowest  vibrational  quantum  and  the 
anharmonicity  constant  equal  to  their  known  spectroscopic  values,  the  depth  of  the  resulting  potential  well  exceeds 
the  experimental  value  of  the  molecule’s  dissociation  energy.  This  fact,  of  course,  means  that  the  Morse  oscillator 
model  is  inaccurate  at  high  vibrational  levels,  but  it  also  gives  an  opportunity  to  formally  model  dissociation  as  a 
bound-bound  transition.  Specifically,  we  will  consider  transitions  to  any  of  the  Morse  oscillator  quantum  levels  -with 
energy  equal  to  or  above  the  experimental  bond  energy  of  the  molecule  to  be  dissociative.  Obviously,  transitions  to 
the  quantum  levels  with  energies  below  the  experimental  bond  energy,  also  approximated  by  the  same  Morse 
oscillator  throughout  the  present  paper,  are  still  interpreted  as  bound-bound  transitions.  The  choice  of  the  ‘cutoff 
bond  energy  is  assumed  to  be  independent  of  the  rotational  energy  of  the  molecules.  The  parameters  of  the 
intermolecular  Morse  potential  used  are  listed  in  Table  I. 


Table  1.  Morse  potential  parameters 


Nitrogen 

Oxygen 

COe,  cm'' 

2359.6 

1580.3 

Xe 

6. 1265  10'^ 

7.639- 10'^ 

D=C0e/4Xe,  cm"' 

96,287 

51,718 

I^exD? 

78,672 

41,268 

Total  number  of  vibrational  levels,  VD=l/2xe 

81 

65 

Number  of  bound  vibrational  levels  (Ev<DexD) 

46 

35 

The  total  probability  of  collisional  dissociation  from  vibrational  level  i  is,  therefore,  a  sum  over  all  final 
“dissociative”  Morse  vibrational  levels: 

(3) 

EjiD 

where  D  is  the  molecule’s  dissociation  (bond)  energy,  and  the  subscript  “VRT'  denotes  vibration-rotation-translation 
energy  transfer. 
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(4) 


Similarly,  for  molecule-molecule  (AB-CD)  collisions: 

AB(i,)+CD{i,)^AB(f,)+CD{f,) 
the  dissociation  probability  can  be  defined  as: 

f\)  (5) 

As  demonstrated  earlier,"’  in  calculating  probabilities  of  vibrational  transitions  in  AB  molecule  at  high 
kinetic  energies  or  temperatures,  vibration-vibration  (VV)  exchange  between  AB  and  CD  can  be  neglected.  Since 
this  work  is  devoted  to  high-temperature  dissociation  processes,  we  will  adopt  this  approximation.  Additionally, 
summing  over  all  possible  vibrational  states  of  CD  molecule  can  be  performed.  The  FHO  transition  probability 
scaling"’  is  then  similar  to  that  in  the  molecule-atom  case: 

'  HT  i' 


all^/,  all;  Q,)  =  i\f\QCUx^{-Q,)Y^ 


■PypT{i-^  A  Q\) 


Here  n  =  min(i,/) ,  and  Qx=AEI  ^  ruj  is  the  dimensionless  energy  (i.e.  the  average  number  of  quanta)  transferred  to 

the  initially  non-vibrating  classical  oscillator  AB  in  an  AB-CD  collision. 

Note  that  in  Refs.  13  and  14,  the  exact  FHO  probabilities  of  Eqs.  (2)  and  (6)  have  been  reduced  to  the 
approximate  analytic  expressions  using  the  asymptotic  behavior  of  the  series  in  the  right-hand-sides  of  Eqs.  (2)  and 
(6)  (e.g.,  see  Ref.  11).  This  permitted  straightforward  analytic  integration  of  the  transition  probabilities  over  the 
collision  energies  and  obtaining  the  closed-form  analytic  expressions  for  the  state-specific  rate  coefficients  for 
vibrational  energy  transfer.  Since  developing  a  fully  analytic  dissociation  model  is  beyond  the  scope  of  the  present 
paper,  we  are  using  the  exact  FHO  scaling  expressions. 


Energy  transfer  to  a  classical  rotating  oscillator  in  3D  collisions;  free-rotation  and  impulsive  approximations 

From  the  previous  section,  it  is  clear  that  to  evaluate  probabilities  of  vibrational  transitions,  including  dissociation, 
an  energy  transferred  to  a  classical  initially  non-vibrating  oscillator  has  to  be  calculated.  For  this,  an  adequate 
potential  energy  surface  should  be  determined,  and  then  classical  equations  of  motion  should  be  integrated  for  the 
three-  or  four-atom  system.  We  will  use  2  approximations  permitting  analytic  solutions. 

III.l  Free-rotation  approximation 


Consider  collisions  between  a  rotating  symmetric  diatomic  molecule  and  an  atom.  For  the  pairwise  atom- 
to-atom  interaction  described  by  repulsive  exponential  functions,  the  atom-molecule  interaction  potential  can  be 
written  as  follows:'^’ 


U (R,  r,  9,  =  2Ae  cosh 


—cos  >9  cos® 

I  2 


(7) 


where  R  is  the  center-of-mass  distance,  r  is  the  separation  of  atoms  in  a  molecule,  9  is  the  rotation  angle  (i.e.  the 
angle  between  R  and  r),  and  <p  is  the  angular  momentum  vector  (i.e.  the  angle  between  the  plane  of  rotation  and  the 
radius  vector  K).  As  in  earlier  work",  we  will  assume  that  the  rotation  is  free,  that  is,  neither  the  magnitude  nor  the 
direction  of  the  angular  momentum  vector  change  in  a  collision.  Then  the  classical  trajectory  can  be  found 
analytically,  and  it  is  essentially  an  exponential  function  of  R(t)  modulated  by  the  periodic  rotational  factor.  The 
resulting  expression  for  the  Q  in  Eq.  (2)  is:" 


Q(E.  s.y.  3. 


te 


sinh" 


nm 


[auy{e,9,^,y)j 


where  E=E^+E,o%  is  the  total  collision  energy,  £=ErJE,  y=AlRj,  b  is  the  impact  parameter,  Rm=2.5  A  is  the  hard- 
sphere  diameter,  5=1 1-/  is  the  number  of  quanta  transferred,  fir=l£(-Ej|/j  h  is  the  average  vibrational  quantum  for 
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the  transition  i^f,  O'  =  ^  ^  ^  q  -  ^  „  _  ^2E/m  ,  m  is  the  collision  reduced  mass,  ttio  is  the  oscillator 

a  k  k 

reduced  mass,  2^=m/mo  is  the  ratio  of  the  collision  and  oscillator  reduced  masses,  and 


r{s,0,(p,y)  =  max 


0, 


(9) 


The  range  of  applicability  of  the  ffee-rotation  approximation  was  discussed  in  the  earlier  paper. In  the 
same  paper,  the  FHO  vibrational  transition  probabilities  calculated  with  this  approximation  were  found  to  agree  very 
well  with  the  full  semiclassical  modeling  using  full  potential  energy  surface.^^ 

An  extension  to  3D  molecule-molecule  collisions  with  rotation  is  straightforward.'''  The  potential  energy 
can  be  written  as: 


U [R,r =  AAe  cosh f — cos i9, cos Icosh f  —  cos 9. 


\  2 


I  2 


COS^ZJj 


(10) 


where  R  is  the  center-of-mass  distance,  r  is  the  separation  of  atoms  in  the  molecules,  i?,  are  the  rotation  angles  (i.e., 
the  angles  between  R  and  each  /•),  and  (pi  are  the  angular  momentum  vectors  (i.e.  the  angles  between  the  planes  of 
rotation  and  the  radius  vector  R).  Having  calculated  the  classical  trajectory  in  this  potential,  assuming  the  free 
rotation,  the  energy  transfer  parameter  Q  turns  out  to  be: ''' 


(  ^>  ^1  >  ^2  >  ^^1  >  *^2  >  ^2  >  T)  ■ 


6. 


cos^  i9j  cos^  q>i^ 


4^ 


sinh' 


TUO, 


;  k=\,  2 


(11) 


Yauy{e„s^,9„(p^,9^,(p^,y)\ 

Here  £=£tr+-£'rot,i+^roi,2  is  the  total  collision  energy,  Ey=Ent,ylE,  y=b^lRj,  b  is  the  impact  parameter,  /?ni=2.5  A  is  the 
hard  sphere  diameter,  4-/kI  is  the  number  of  quanta  transferred,  o\=  1  E^yE^x  I  ^  is  the  average  vibrational 


quantum  for  the  transition  4— ^/k,  = 


■  ha>.  /.  _ ,  .  ,  , 

- - — ,  0^  =  — -  ,  M  =  v2£ /  w  ,  m  IS  the  collision  reduced  mass,  and 


a  k 


r{e„£^,9„(p^,9^,(p^,y)  =  max 


^,__sin2a  cost,,  ^_sm2.9  .osy, 


(12) 


III.2  Impulsive  (sudden)  approximation 


In  this  approximation,  AB  is  modeled  by  two  hard  spheres,  A  and  B,  connected  by  an  initially  non- 
deformed  spring.  Collision  with  an  atom  C,  also  modeled  as  a  hard  sphere,  reduces  to  a  hard-sphere  collision  of  B 
and  C.  The  center-of-mass  kinetic  energy,  AB  rotational  energy,  and  the  set  of  orientation  angles  at  the  moment  of 
collision  fully  define  the  outcome  of  the  collision.  Applying  energy  and  momentum  conservation,  the  post-collision 
component  of  the  velocity  of  B  atom  along  the  A-B  axis  can  be  determined,  and  the  post-collision  vibrational  energy 
of  the  AB  molecule  calculated.  Similarly,  AB-CD  collisions  can  be  analyzed  in  the  assumption  that  only  hard 
spheres  B  and  C  collide  directly. 

The  hard-sphere  collision  of  B  and  C  atoms  can  be  treated  as  being  either  “smooth-sphere”,  where  the 
collision  is  elastic  along  the  line  connecting  centers  of  B  and  C,  with  no  forces  normal  to  that  line,  or  “rough-sphere” 
one,  where  the  collision  is  elastic  along  the  line  coimecting  B  and  C  centers  and  fully  inelastic  perpendicular  to  that 
line. 

When,  in  a  collision,  B  atom  receives  a  velocity  increment  Av  along  the  A-B  axis,  the  initially  non¬ 
vibrating  AB  molecule  acquires  vibrational  energy 

where  m  is  the  mass  of  each  of  the  atoms  A  and  B.  Since,  in  the  present  model,  B  collides  directly  only  with  C,  it  is 
convenient  to  introduce  B  velocity  increments  along  the  B-C  axis,  Avp ,  and  perpendicular  to  the  B-C  axis  and 

parallel  to  the  A-B-C  plane,  Av^^ ,  so  that 

Av  =  Avp  cos  0  +  AVj^  sin  ^  .  (14) 

Here  6  is  the  angle  between  the  A-B  molecular  axis  and  the  B-C  axis;  ^  e  [O,  ;r/2] .  Note  that  the  B  atom  velocity 
increment  component  perpendicular  to  both  B-C  axis  and  A-B-C  plane  does  not  contribute  to  Av . 
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Consider  now  an  AB-C  collision  in  the  AB-C  center-of-mass  system.  Translational  velocities  of  the 
colliding  atoms  B  and  C,  and  ,  antiparallel  to  each  other,  are  related  to  the  center-of-mass  kinetic  energy, 
E„  of  the  colliding  pair  by  the  momentum  and  energy  equations: 

1  2  1  2 

+-wVc,„  =£, ,  (15) 

where  M  is  the  mass  of  C  atom.  The  direction  of  and  ,,,  can  be  characterized  by  a  ,  the  angle  between  the 
velocity  of  the  C  atom  (in  the  AB-C  center-of-mass  system)  and  the  B-C  axis,  and  by  p ,  the  angle  between 
^a./r  or  ^c.B-  ®od  the  straight  line  parallel  to  the  A-B-C  plane  and  normal  to  the  B-C  axis.  The  intervals  of  variation 
of  these  angles  are:  ar€[0,/r/2]  and  <p€.\Kl2-a,  7r-a\. 

Additionally,  atom  B  has  a  rotational  velocity,  v^  ,  whose  magnitude  can  be  expressed  in  terms  of  the 
initial  rotational  energy  of  AB,  £): 


2~mv 


B,rot  5 


and  whose  direction  (obviously,  normal  to  the  A-B  axis)  can  be  characterized  by  the  angle  between  v^  and  the 
B-C  axis;  ^g[^/2-0,  nl2^B\. 

Considering  an  elastic  collision  of  B  and  C  atoms  along  the  B-C  line,  from  the  energy  and  momentum 
conservation  equations  we  obtain: 

*  2M  r/  ^  2M  \  \(M  +  2m)E,  [e  J 

^v_= -  (v„  „ +Vr  ,  lcOSa  +  V„  COS;0  = -  - ^—^cosa  +  J—cos0  (17) 

Here  we  expressed  velocities  v^  ,  Vg  and  v^  „  through  the  energies  Er  and  E,  according  to  equations  (15)  and 
(16).  Since  in  the  smooth  spheres  case  AVj^  =  0 ,  substituting  Eqn.  (17)  into  Eqn.  (14),  and  inserting  the  result  into 
Eqn.  (13),  we  obtain  the  energy  transferred  to  the  initially  non-vibrating  oscillator  in  the  smooth  spheres,  molecule- 
atom  case: 


|(yt?^)Acos«  +  V^cos,glcosg] 

-f  +  m  V  M  V  r  y 

\  y 


In  the  rough  spheres  case,  one  needs  to  calculate  projections  of  all  velocities  of  B  and  C  atoms  onto  an  axis 
perpendicular  to  the  B-C  line  and  parallel  to  the  A-B-C  plane.  Then,  from  the  energy  and  momentum  equations  for 
an  absolutely  inelastic  collision  of  B  and  C  atoms  along  that  axis,  one  obtains: 

Mr/  \  „  ^ _ 


Av,  =- 


M  +  m 


( „  +  Vc,  „  )  cos  97  -  v^  cos  cos  61  /  sin  0 


Expressing  velocities  v,  ,  Vg  andv^  ,,,  through  the  energies  Er  and  E,  according  to  equations  (15)  and  (16), 

then  substituting  Eqns.  (17)  and  (19)  into  Eqn.  (14),  and  inserting  the  resulting  formula  for  Av  into  Eqn.  (13),  we 
obtain  the  energy  transferred  to  the  initially  non-vibrating  oscillator  in  the  rough  spheres,  molecule-atom  case: 

(2cosgcos^H- cos sin +  ^Je^ cos ^ cos 0  |  .  (20) 

Similar  derivations  can  be  performed  for  molecule-molecule  collisions.  When  writing  energy  and  momentum 
conservation  equations  for  B-C  collisions,  additional  (compared  with  the  molecule-atom  case)  C  atom  velocity 
components,  corresponding  to  rotational  and  vibrational  motion  of  CD  molecule,  have  to  be  included,  leading  to 
additional  terms  in  equations. 

In  the  molecule-molecule,  smooth  spheres  case: 


M  +  m 


M  I  I(M  +  m)E, 


cosa  + 


^cos/3- 


E^  cos  T]  cos  /,+J — Er  cos  cos  0 


Here,  of  course,  we  consider  AB-CD  collisions,  with  B  and  C  atoms  actually  colliding.  The  mass  of  each  of  the 
atoms  C  and  D  equals  M.  E^  and  E\.  are  the  pre-collision  vibrational  and  rotational  energies  of  the  CD  molecule. 


99 


Arguments  of  the  new  cosines  are;  rj  -  the  phase  of  CD  oscillations  in  the  moment  of  collision,  -  the  angle 
between  B-C  and  C-D  axes,  and  -  the  angle  between  the  rotational  velocity  of  C  atom  and  the  B-C  axis.  The 

intervals  of  variation  of  the  new  arguments  are; 

;7e[0,  27t\,  e[0, 

^2  e  [n- /  2  -  y, ,  ;t/2  +  y,  ],  if  y,  e  [O,  nl2\,  (22) 

e[y, -;t/2,  3;r/2-y,],  ify,  e(;r/2,  n\ 

In  the  molecule-molecule,  rough  spheres  case; 
f 

1 


AE: 


M 


2  M  +  m 


(M  +  m)E,  .  ^ 

—  (2  cos  a  cosy  +  cos  ^ sin  0) 


3  '  P*ii"  I 

cosj0cos0+y— cos77(2cosy|Cos^  +  cos^,  sin 


(23) 


+  Er{^  COS  Yi  cos  9  +  cos  sin  6) 

The  new  angles,  (p^  and  <p^ ,  are  the  angles  between  vectors  of  vibrational  and  rotational  velocity,  respectively,  of  C 
atom  and  the  straight  line  parallel  to  the  A-B-C  plane  and  perpendicular  to  the  B-C  axis.  The  intervals  of  their 
variation  are; 

((\e{nl2-y^,  7rl2  +  y^],  ify|€[0,  ;r/2], 
e(y,  -nl2,  3;T/2-y,],  if  y,  &{nl2, 
and 

(p^B[7:l2-y^,  7cl2  +  y^],  ify2€[0,  ;r/2], 

&[y^-nl2,  2nl2-y^\,  if y^  e[7tl2,  n\ 


(24) 


Dissociation  probabilities,  cross  sections,  and  rate  coefficients 

Formulas  (8),  (9),  (1 1),  (12),  and  (18)-(24)  for  Q  and  AE  ,  when  inserted  into  FHO  scaling  formulas  (2)  and  (6), 
give  sets  of  state-to-state  transfer  probabilities.  Summing  over  the  final  states  according  to  Eqs.  (3)  and  (5), 
dissociation  probabilities  are  found  as  functions  of  vibrational,  rotational,  and  translational  energies,  and  of  the  set  of 
angles.  Integrating  over  the  angles,  and  assuming  rotational-translational  equilibrium  with  temperature  T,  gives  the 
dissociation  cross  section  as  a  function  of  vibrational  quantum  number  i,  total  (translational  plus  rotational)  collision 
energy  E,  and  the  temperature  T. 

The  dissociation  cross-section  for  molecule-atom  collisions  in  the  free-rotation  approximation  is; 

iml  111  K  ^  0  \j-j\  0  0  /  ^25) 

\\dy\dJ\d9^d(p  Y,Pij{E,  £,  9,  (p,  y) 

V-'  /  ^  0  0  0  0  f 

where  1  is  the  moment  of  inertia  of  the  AB  molecule.  For  molecule-molecule  collisions  in  the  free-rotation 
approximation; 
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<^a,Ai^E,T)  = 


rax  vlmx  /2nax 

|rf/(2/  +  l)  I  4(2/, +1)  J  4(2/, +1) 


K  n  K  n 

■  _  d&^  IdVi  Jrf«9j  ■Y.Pif 


7  —  J4  p>9, 1^9’. 


/ 

where  /,  and  /^  are  the  moments  of  inertia  of  AB  and  CD  molecules,  respectively.  Similar  formulas  can  be  written 
for  the  impulsive  approximation.  In  Eqs.  (25)  and  (26),  the  rigid  rotor  approximation  has  been  used.  The  limits  of 
integration  over  the  angular  momenta  in  Eqs.  (25)  and  (26)  are  determined  from  the  condition  that  the  cross  section 
is  evaluated  at  the  constant  total  energy  E  of  the  classical  degrees  of  freedom  (i.e.  the  sum  of  translational  energies 
of  the  radial  and  the  orbital  motion  and  the  rotational  energies,  E=Erad+Eorb+Eroi,k).^'’  The  approximate  integration 
forms  in  the  right-hand  sides  of  Eqs.  (25)  and  (26)  are  obtained  by  replacing  integration  over  the  rotational  and 
orbital  angular  momenta  by  integration  over  the  dimensionless  rotational  energies,  E^Ef^  ^JE-,  and  the  dimensionless 
impact  parameter  y=b^lRj'.  In  this  procedure,  replacing  the  appropriate  integration  linuts  for  and  y  by  unity  does 
not  significantly  affect  the  integration  accuracy  since  the  dissociation  cross  section,  oa,ss,  was  found  to  steeply  drop 
both  at  the  high  values  of  rotational  energies  and  the  impact  parameter  (i.e.  at  the  low  values  of  the  translational 
energy  of  the  radial  motion). 

In  this  work,  cross-sections  were  calculated  at  60  values  of  total  collision  energy  in  the  range  10^-10®  cm'*. 
For  each  total  collision  energy  E,  1000  randomly  picked  phase  space  points  were  used  for  evaluation  of  multi¬ 
dimensional  cross-section  integral.  Previous  calculations  using  this  approach  by  Billing^'*’  showed  that  the  use  of 
only  100-400  phase  space  points  is  sufficient  to  obtain  a  20-30%  confidence  in  the  averaged  cross  sections.  In  the 
present  calculations,  the  use  of  1000  phase  space  points  allows  about  10-20%  accuracy  in  the  cross  sections.  To 
account  for  the  oscillator  reaction  on  the  classical  trajectory,  or,  in  other  words,  to  enforce  the  detailed  balance,  a 
symmetrization  procedure^'’  was  used,  whereby  the  effective  energy  U  was  substituted  for  E,  corresponding  to  the 
arithmetically  averaging  velocity  before  and  after  collision: 

A£,,  (AE,.)' 

U  =  E  +  — E  +  \ — ’U-.  £■  =A£:,.,/4  (27) 

2  16£  ^  ’ 

Averaging  the  dissociation  flux  over  Maxwellian  energy  distribution  gives  vibrational  state-specific 
dissociation  rate  coefficient  (v  is  the  vibrational  quantum  number): 

J"  K.(vV,r)exp[-^).(^]  (28, 

The  integral  in  Eq.  (28)  is  evaluated  numerically  by  the  RATECONS  procedure*'  using  all  60  values  of  the 
cross  sections  given  by  Eqs.  (25)  and  (26).  The  estimated  accuracy  of  numerical  integration  is  about  10%.  Further 
averaging  of  the  state-specific  rate  coefficients  over  the  vibrational  energy  mode  distribution  function  would  give 
the  total  dissociation  rate  coefficient.  In  the  present  paper,  this  is  done  for  both  nonequilibrium  and  equilibrium 
vibrational  mode  energy  distributions.  In  both  cases,  the  energy  distributions  are  assumed  to  be  Boltzmaiui  with 
vibrational  temperature  Ty  and  the  gas  temperature  T,  respectively.  In  the  latter  case,  the  thermal  dissociation  rate 
constant  is  obtained. 


E-nm  ~  EE  if  /  4 


4„(v.T)=  -  J«7<„>,(v,t/,T)exp  - 


in.  Results  and  discussion 

To  check  the  difference  between  vibrational  energy  transfer  predicted  by  free-rotation  and  impulsive 
approximations,  the  parameter  Q  =  AE  /  ha) ,  averaged  over  collision  angles,  was  calculated  for  02-Ar  and  Oa-O, 
collisions  and  plotted  in  Fig.  1  (a,  b)  as  a  function  of  the  total  (translational  plus  rotational)  collision  energy.  In  all 
free-rotation  calculations  done  in  this  work,  the  parameters  of  the  exponential  repulsive  potential  of  Eq.  (7)  were 
taken  to  be  a  =4.0x10*  cm"' ,  which  is  the  commonly  accepted  number  for  nitrogen  and  oxygen,'®''*'  **’  and 
A=1730  eV.*''  As  seen  in  Fig.  1  (a,  b),  there  is  virtually  no  difference  between  smooth  and  rough  sphere  impulsive 
approximations.  At  very  high  energies,  free-rotation  and  impulsive  approximations  converge  to  close  asymptotics, 
which  is  physically  reasonable  and  provides  a  consistency  check.  At  moderate  and  low  energies,  the  impulsive 
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approximation  is  seen  to  overpredict  the  energy  transfer  by  orders  of  magnitude.  Again,  this  could  be  expected, 
since  in  the  impulsive  approximation  collision  duration  is  assumed  to  be  much  shorter  than  the  vibrational  period, 
while  low-velocity  collisions  last  much  longer  than  the  period  of  vibrations,  thus  reducing  the  efficiency  of  energy 
transfer. 

Fig.  2  (a,  b)  shows  the  dissociation  probability  from  ground  vibrational  state  versus  total  collision  energy. 
Similar  to  Fig.  1,  free-rotation  and  impulsive  results  are  close  to  each  other  at  high  energies,  while  at  low  energies 
impulsive  approximation  predicts  an  excessively  high  dissociation  probability. 

An  important  step  in  validating  any  dissociation  model  is  to  compare  predicted  dissociation  rates  in  thermal 
equilibrium  with  experimental  data.  Strictly  speaking,  this  would  require  full  master  equation  modeling,  with 
dissociation,  recombination,  and  energy  transfer  processes  all  taken  into  account.  Such  modeling  can  be  done  using 
the  newly  derived  dissociation  rates  and  FHO-based  vibrational  energy  transfer  rates,’®’  ''''  similar  to  the 
modeling  performed  earlier.”’  However,  such  modeling  is  beyond  the  scope  of  this  paper.  For  the  first-cut 
comparison,  in  this  work  we  calculated  thermal  dissociation  rates  for  nitrogen  and  oxygen  assuming  the  equilibrium 
Boltzmaiui  vibrational  energy  distribution  function.  Of  course,  if  atom-atom  recombination  is  weak,  the  irreversible 
dissociation  would  depopulate  very  high  vibrational  states,  reducing  the  dissociation  rate  by  perhaps  a  factor  of  2- 
3  1’  7, 33  ^  other  hand,  recombination  is  significant,  highly  excited  states  would  be  close  to  equilibrium.  Thus, 

a  caution  is  needed  when  comparing  computed  thermal  dissociation  rates  with  experiments,  as  the  degree  to  which 
highly  excited  states  were  close  to  thermal  equilibrium  in  a  specific  experiment  is  often  not  clear. 

In  Fig.  3,  a  comparison  with  the  dissociation  rate  inferred  from  shock  tube  measurements  is  shown  for  O2- 
Ar  collisions.  Both  free-rotation  and  impulsive  energy  transfer  models  in  conjunction  with  FHO  scaling  give  an 
excellent  agreement  with  experiments  in  a  very  wide  temperature  range,  especially  considering  that  the  experimental 
fit  has  an  uncertainty  range  within  a  factor  of  2-3.  The  agreement  of  the  impulsive  model  with  both  experimental 
data  and  free-rotation  predictions  may  seem  somewhat  surprising.  The  explanation  lies  in  the  fact  that  at  thermal 
equilibrium  the  dissociation  proceeds  chiefly  from  very  high  vibrational  levels,  close  to  the  dissociation  threshold. 
This  is  well  known  from  the  literature,’  and  is  confirmed  in  our  calculations  (see  below).  Since  spacing  between 
vibrational  levels  is  very  small  near  the  dissociation  threshold,  the  Massey  parameter,’  that  is,  the  product  of 
vibrational  frequency  and  the  collision  duration,  is  small  even  for  particles  moving  with  mean  thermal  velocity.  This 
means  that  most  dissociative  collisions  at  thermal  equilibrium  can  be  adequately  described  in  the  impulsive 
approximation. 

Figure  4  shows  calculated  and  experimental  thermal  dissociation  rate  coefficients  in  oxygen.  One  can  see 
that  both  models  somewhat  underpredict  the  dissociation  rate  compared  to  the  experimental  data, 

2.2xl0’*exp(-59380/r)[l-exp(-2240/r)]  crn^  I  [mole  ■  s) .  The  impulsive  model,  which  provides  for  more 

efficient  rotational  energy  transfer  into  the  dissociating  bond  than  the  free-rotation  model  does,  appears  to  be  closer 
to  experiments  in  this  case.  It  is  also  instmctive  to  consider  the  ratio  of  thermal  dissociation  rate  coefficients  in  O2- 
O2  and  02-Ar  collisions.  Both  free-rotation  and  impulsive  models  predict  the  ratio  to  be  about  3  and  to  depend  on 
temperature  very  weakly.  That  molecule-molecule  collisions  are  more  effective  than  molecule-atom  ones 
demonstrates  the  role  of  rotational  energy  in  dissociation  and  is  qualitatively  consistent  with  theoretical 
expectations.’  Nikitin’  argued  that  the  vibrational  energy  E*,  that  approximately  separates  the  regions  of  single¬ 
quantum  and  multiquantum  (including  bound-free)  transitions  is  determined  by  the  criterion  that  the  Massey 
parameter  calculated  with  vibrational  frequency  at  the  energy  level  be  of  the  order  of  unity.  On  the  basis  of  this 
criterion,  Nikitin’  estimated  reduction  in  E*  due  to  efficient  transfer  of  rotational  energy  of  the  non-dissociating 
molecule  into  vibrational  mode  of  the  highly-excited  dissociating  molecule  and  predicted  the  ratio  of  thermal 
dissociation  rate  coefficients  in  O2-O2  and  02-Ar  collisions  to  be  about  10-20,  independent  of  temperature. 
Experimental  data  quoted  by  Stupochenko,  Losev,  and  Osipov,^^  specifically,  O2-O2  data^'*'^®  and  02-Ar  data^’  seem 
to  indicate  that  the  dissociation  rate  ratio  decreases  from  30-40  at  3500  K  to  about  5-10  at  7000  K.  Since  at  thermal 
equilibrium  dissociation  in  both  O2-O2  and  02-Ar  collisions  is  determined  by  the  processes  at  very  high  vibrational 
states,  one  possible  reason  for  the  less  than  perfect  agreement  between  our  predictions  and  the  experiment  could  be 
an  error  introduced  by  modeling  the  highly  excited  states  by  Morse  oscillator.  Further  refinement  of  the  theoretical 
model,  especially  with  regard  to  the  energy  spectmm  at  high  vibrational  levels,  would  clarify  the  issue.  However, 
experimental  data  of  different  authors  also  have  to  be  treated  with  caution.  For  example,  Losev  et.  al.*’  based  on 
experimental  data  of  Losev  and  Shatalov,^”  recommend  the  dissociation  rate  value  of 
2.2x10'' exp(-59380/7)[l  -exp (-2240/7’)]  cm^ l{mole-s^  for  O2-O2  collisions.  Together  with  the  02-Ar  rate 

6.0x10”  (59380/r)exp(-59380/r)  err?  l{mole-s)  of  Camac  and  Vaughan,^’  this  would  give  the  ratio  of  thermal 
dissociation  rate  coefficients  in  O2-O2  and  02-Ar  collisions  of  about  10,  weakly  dependent  on  temperature  (see  Fig. 
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4).  Considering  the  fact  that  both  rates  are  known  within  a  factor  of  3  at  best,  this  would  mean  a  reasonably  good 
agreement  between  our  model  and  experiments.  A  detailed  analysis  of  the  experiments  quoted,  as  well  as  further 
shock-tube  experimental  studies,  could  be  helpful  in  developing  a  better  understanding  of  dissociation  processes. 

Figures  5  and  6  depict  calculated  and  experimental  thermal  dissociation  rate  coefficients  for  Na-Ar  and  N2- 
N2  collisions.  As  seen  from  the  figures,  the  agreement,  within  a  factor  of  2-3  experimental  uncertainty,  is  uniformly 
excellent  over  the  wide  temperature  range.  Interestingly,  the  calculated  ratio  of  the  thermal  dissociation  rate 
coefficients  for  N2-N2  and  N2-Ar  collisions  is  about  3,  weakly  changing  with  temperature,  which  is  in  a  very  good 
agreement  with  experiments. 

As  is  well  known,’’  rotational  relaxation  behind  shocks  is  quite  fast,  while  vibrational  relaxation  is 
much  slower.  Thus,  dissociation  rate  coefficient  as  a  function  of  vibrational  energy  and  translational-rotational 
temperature  is  quite  useful  for  shock  modeling.  When  averaged  over  nonequilibrium  vibrational  distribution,  such 
vibrational  energy-dependent  rate  coefficient  could  be  converted  into  a  two-temperature  rate  coefficient  k{T^,T) . 
The  calculated  vibrational  energy-dependent  rates  normalized  by  the  dissociation  rate  from  the  uppermost 
vibrational  level  at  several  temperatures  are  shown  in  Fig.  7  for  N2-Ar  collisions.  Very  steep  exponential  increase  of 
the  rates  with  vibrational  energy  is  clearly  seen.  Perhaps  more  informative  is  Fig.  8,  where  the  rates  from  Fig.  7  have 
been  multiplied  by  the  equilibrium  Boltzmaim  factor  =  exp(-£„  /  T) .  The  plotted  quantity,  therefore,  represents 

contributions  of  various  vibrational  states  to  the  total  thermal  dissociation  rate.  As  expected  from  theory,’  the 
principal  contribution  to  thermal  dissociation  comes  from  states  with  vibrational  energy  close  to  the  dissociation 
threshold. 

The  plot  in  Fig.  8  is  helpful  in  understanding  the  degree  of  preferentiality  of  dissociation.  The  well-known 
Marrone-Treanor  model’  assumes  that  the  probability  of  dissociation  from  vibrational  state  with  energy  Ey  becomes 
non-zero  when  translational  collision  energy  exceeds  the  threshold  (£>-£,),  and  that  there  is  an  additional 


exponential  increase  of  the  probability  with  vibrational  energy,  so  that  the  vibrational  energy-dependent  rate  is 
proportional  to 


*wr(^v.7’)«exp  — 


- =7^ 

•exp 

1  T  J 

(29) 


where  the  empirical  parameter  Ui^n-  is  a  measure  of  preferentiality  of  the  dissociation  from  higher  vibrational  states. 
The  smaller  Umt,  the  more  preferential  dissociation  is,  while  =  00  conesponds  to  non-preferential  dissociation. 


Multiplying  the  Marrone-Treanor  rate  [Eq.  (29)]  by  the  Boltzmann  factor  exp(-£,/7’) ,  normalizing,  and  plotting 


the  product  in  a  semi-logarithmic  scale  would  give  a  linear  function  of  vibrational  energy,  with  the  slope 
determining  the  preferentiality  parameter.  In  Fig.  8  that  shows  such  a  plot,  the  free-rotation  approximation  is  seen  to 
predict  dissociation  rates  that  are  very  preferential  towards  higher  vibrational  states.  However,  interpreting  the  rates 
in  terms  of  a  single  constant  parameter  Umt  is  quite  difficult.  First,  the  functions  plotted  in  a  senh-logarithmic  scale 
in  Fig.  8  deviate  substantially  from  linearity.  Second,  slopes  of  the  curves  at  each  vibrational  energy  seem  to 
decrease  with  increasing  temperature,  so  that  the  effective  parameter  Umt  would  increase  with  temperature,  making 
dissociation  less  preferential  at  higher  temperatures.  The  impulsive  model,  as  seen  in  Fig.  8,  gives  rates  that  are  less 
preferential  than  those  given  by  the  free-rotation  model,  so  that  at  very  high  temperatures  the  impulsive-model  rates 
are  virtually  non-preferential. 

Finally,  the  new  model  was  tested  against  experimental  data  on  nonequilibrium  dissociation  of  nitrogen.  Losev 
et.  al*’  *’  have  measured  both  vibrational  and  translational-rotational  temperatures,  and  the  nonequilibrium 
factor  k(Tv,T)/k(T),  i.e.  the  ratio  of  nonequilibrium  and  thermally-equilibrium  dissociation  rates,  has  been  inferred 
from  shock-tube  experiments.  In  Ref.  38,  the  best  fit  to  the  measured  nonequilibrium  factor  is  obtained  using  a  two- 
temperature  dissociation  model  *, 

A(r^,r)_i-exp(-g/r,) 
kiT)  l-exp(-6'/r)  ^ 

with  P=3.  In  Eq.  (30),  0=3353  K  is  the  characteristic  vibrational  temperature  of  N2.  As  shown  in  Fig.  9,  the  present 
model  agrees  well  with  the  data  inferred  from  the  experiments  (Eqn.  (30)  vrith  P=3),  considering  that  the  latter  are 
known  within  a  factor  of  2-3.  The  theoretical  two-temperature  dissociation  rate  k(Tv,T)  shown  in  Fig.  9  is  obtained 
by  averaging  the  state-specific  rates  k(Ev,T)  over  the  nonequilibrium  Boltzmann  distribution  over  the  vibrational 
mode  energy,  with  the  vibrational  temperature  Ty.  Note  that  Figs.  8  and  9  plot  two  essentially  different 
characteristics  of  the  dissociation  process.  Particularly,  Fig.  8  shows  the  relative  contribution  of  each  vibrational 
level  into  the  overall  dissociation  rate  at  thermal  equilibrium,  k(Ev,T)feq(Ev,T)/k(T),  where  fe,(Ey,T)  is  the 
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equilibrium  relative  population.  On  the  other  hand.  Fig.  9  shows  the  sum  of  contributions  of  all  vibrational  levels 
into  the  overall  dissociation  rate  at  nonequlibrium  conditions,  at  Tv<T.  In  fact,  only  the  asymptotic  limit  of  the  ratio 
k(Tv,T)/k(T)  at  Tv->0,  plotted  in  Fig.  9  shows  the  relative  contribution  of  vibrational  level  v=0  alone,  i.e. 
k(Ev=0,T)/k(T). 

IV.  Conclusions 

The  new  semiclassical  theory  of  dissociation  of  diatomic  molecules,  suggested  in  this  paper,  is  based  upon  an 
anharmonicity-corrected  and  energy-symmetrized  Forced  Harmonic  Oscillator  (FHO)  quantum  scaling  in 
conjunction  with  free-rotation  or  impulsive  energy-transfer  models.  Although  the  theory  is  not  fully  analytical,  it  is 
computationally  simple,  and  it  can  predict  both  state-specific  and  thermal  dissociation  rates,  explicitly  accounting 
for  molecular  rotation  and  three-dimensional  collisions,  without  any  adjustable  parameters. 

Calculated  thermal  dissociation  rates  of  nitrogen  and  oxygen  are  in  excellent  agreement  with  shock-tube  data 
over  a  wide  temperature  range.  Agreement  with  experimental  data  on  nonequilibrium  dissociation  of  nitrogen  is  also 
quite  good. 

Predicted  vibrational  energy-dependent  rates  show  that  dissociation  from  very  high  vibrational  states  is  strongly 
favored  in  the  free-rotation  approximation.  However,  the  preferential  nature  of  dissociation  is  difficult  to  interpret  in 
terms  of  a  single  Marrone-Treanor  parameter  Umt,  as  this  parameter  would  be  vibrational  energy  dependent  and 
increase  with  temperature.  The  impulsive  approximation  for  energy  transfer  results  in  dissociation  rates  with 
relatively  little  preference  for  high  vibrational  states. 

A  possible  future  refinement  of  the  theory  would  involve  using  a  more  sophisticated  and  precise  energy 
spectmm  of  highly  vibrationally  excited  molecules.  For  practical  application  in  CFD  or  DSMC  analyses,  the  theory 
could  be  used  to  either  generate  a  set  of  tabulated  rates  or  to  develop  simple  interpolation  formulas. 


REFERENCES 

’  C.Park,  Nonequilibrium  Hypersonic  Aerothermodynamics  (Wiley,  N.Y.,  1990). 

^  S.O.Macheret,  Paper  AIAA  97-2501,  1997. 

^  P.V.Marrone  and  C.E.Treanor,  Phys.  Fluids  6,  1215  (1963). 

“*  S.A. Losev,  Paper  AIAA  96-2026, 1996. 

^  A.L.Sergievskaya,  E.A.Kovach,  S.A.Losev,  and  N.M.Kuznetsov,  Paper  AIAA  96-1895,  1996. 

*  C.Park,  J.  Thermophys.  and  Heat  Transfer  2,  8  (1988);  ibid.,  3,  233  (1989). 

^  E.E.Nikitin,  Theory  of  Elementary  Atomic  and  Molecular  Processes  in  Gases  (Clarendon  Press,  Oxford,  1974), 
Chap.  7. 

*R.N. Schwartz,  Z.I.Slawsky,  and  K.F.Herzfeld,  J.  Chem.  Phys.  22,  767  (1954). 

®S.P.Sharma,  W.Huo,  and  C.Park,  Paper  AIAA  88-2714,  1988. 

I.V.Adamovich,  S.O.Macheret,  J.W.Rich,  and  C.E.Treanor,  AIAA  Journal  33,  1064  (1995). 

”  I.V.Adamovich,  S.O.Macheret,  C.E.Treanor,  and  J.W.Rich,  J.  Thermophys.  and,Heat  Transfer  12,  57  (1998). 
LV.Adamovich,  J.W.Rich,  and  S.O.Macheret,  J.  Thermophys.  and  Heat  Transfer  11,  261  (1997). 

I. V.Adamovich  and  J.W.Rich,  J.  Chem.  Phys.  109,  7711  (1998). 

'''I.V.Adamovich  and  J.W.Rich,  Paper  AIAA  99-3565, 1999 

S.O.Macheret  and  J.W.Rich,  Chem.  Phys.  174,  25  (1993). 

a)  S.O.Macheret,  A.A.Fridman,  I.V.Adamovich,  J.W.Rich,  and  C.E.Treanor,  Paper  AIAA  94-1984,  1994;.b) 
I.V.Adamovich,  S.O.Macheret,  J.W.Rich,  C.E.Treanor,  and  A.A.Fridman,  “Vibrational  Relaxation,  Nonequilibrium 
Chemical  Reactions,  and  Kinetics  of  NO  Formation  Behind  Strong  Shock  Waves”,  in  M.  Capitelli  (ed.)  “Molecular 
Physics  and  Hypersonic  Flows”,  NATO  Advanced  Study  Institute  Series,  vol.  482,  Kluwer,  1996,  p.  85-104 
E.H.Kemer,  Can.  J.  Phys.  36,  371  (1958). 

'*  C.E.Treanor,  J.  Chem.  Phys.  43,  532  (1965). 

A.Zelechow,  D.Rapp,  and  T.E.Sharp,  J.  Chem.  Phys.  49,  286  (1968). 

J. D.Kelley,  J.  Chem.  Phys.  56,  6108  (1972). 

A.J.Banks,  D.C. Clary,  and  H.-J.Werner,  J.  Chem.  Phys.  84,  3788  (1986). 

^^G.D.Billing,  Chem.  Phys.  107, 39  (1986). 

M.M.Maricq,  E. A. Gregory,  C.T.Wickham-Jones,  D.J.Cartwright,  and  C.J.S.M. Simpson,  Chem.  Phys.  75,  347 
(1983). 

G.D.Billing  and  E.R.Fisher,  Chem.  Phys.  43,  395  (1979). 


104 


G.D.Billing  and  R.E.Kolesnick,  Chem.  Phys.  Lett.  200, 382  (1992). 

R. Hemandez,  R.Toumi,  and  D.C.CIary,  J.  Chem.  Phys.  102,  9544  (1995). 

J.M.Price,  J.A.Mack,  C.A.Rogaski,  and  A.M.Wodtke,  Chem.  Phys.  175,  83  (1993) 

^*H.Park  and  T.G.Slanger,  J.  Chem.  Phys.  100, 287  (1994). 

"M.Klatt,  l.W.M.Smith,  R.P.Tuckett,  and  G.N.Ward,  Chem.  Phys.  Lett.  224, 253  (1994) 

^“G.D.Billing,  in  Nonequilibrium  Vibrational  Kinetics,  ed.  by  M.Capitelli  (Springer,  Berlin,  1986),  p.  85. 
G.D.BiIling,  Computer  Physics  Communications  32,  45  (1984) 

G.D.BilIing,  Computer  Physics  Communications  44, 121  (1987) 

Ye.V.Stupochenko,  S.A.Losev,  and  A.I.Osipov,  Relaxation  in  Shock  Waves  (Springer-Verlag,  New  York,  1967). 

”  S.A.Losev,  Doklady  Acad.  Sci.  SSSR  141,  894  (1961). 

”  D.L.Matthews,  Phys.  Fluids  2,  1 70  (1959). 

“S.R.Byron,  J.  Chem.  Phys.  30,  1380  (1959). 

”M.Camac  and  A.Vaughan,  J.  Chem.  Phys.34,  460  (1961). 

S. A.Losev  and  G.G.Chemyi,  eds.,  “Analysis  of  Russian  Experimental  Investigations  for  Validation  of  Kinetic 
Models  for  Hypersonic  Flows  in  Thermal  Nonequilibrium,”  Final  Report  on  the  Research  Grant  No.  93-0987-01  by 
and  between  the  North  Carolina  State  University  and  Avogadro  Center  of  Institute  of  Mechanics  (Moscow  State 
University,  1993). 

^“S.A.Losev  and  O.P.Shatalov,  Doklady  Acad.  Nauk  SSSR  185, 293  (1969)  [Russian]. 

J.P. Appleton,  M.Steinberg,  and  D.J.Liquomik,  J.  Chem.  Phys.  48,  599  (1968). 

S.Byron,  J.  Chem.  Phys.  44,  1378  (1966). 

S.A.Losev  and  M.S.Yalovik,  Khimiia  Vysokikh  Energii  (High  Energy  Chemistry)  4, 202  (1970)  [Russian]. 


Fig. la,  b.  Average  number  of  quanta  transferred  to  initially  non-vibrating  oxygen  molecule  versus  total  collision 
energy,  as  given  by  impulsive  and  free-rotation  models,  (a)  02-Ar  collisions;  (b)  O2-O2  collisions 
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Fig.  2  a,  b.  Probability  of  oxygen  molecule  dissociation  from  the  ground  vibrational  state  versus  total  collision 
energy,  as  calculated  with  impulsive  and  free-rotation  models,  (a)  02-Ar  collisions;  (b)  O2-O2  collisions 
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Fig.  3.  Computed  and  inferred  from  experiments^’  thermal  dissociation  rate  constant  for  oxygen  diluted  iii  argon 
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Fig.  4.  Computed  and  inferred  from  experiments^'*"’*  thermal  dissociation  rate  constant  for  oxygen 
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Fig.  5.  Calculated  and  experimental''*’  nitrogen  thermal  dissociation  rates  in  N2-Ar  collisions  versus  temperature. 


k(T),  cm^/sec 


Temperature,  K 

Fig.  6.  Calculated  and  experimental'"  thermal  dissociation  rate  of  nitrogen  in  N2-N2  collisions  versus  temperature. 
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Fig.  7.  Calculated  state-specific  dissociation  rates  in  N2-Ar  collisions,  normalized  by  the  dissociation  rate  from  the 
uppermost  vibrational  level. 
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Fig.  8.  Calculated  state-specific  dissociation  rates  in  N2-Ar  collisions,  multiplied  by  the  Boltzmann  population  factor 
and  normalized  by  the  thermal  dissociation  rate. 
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Fig.  9.  Comparison  of  calculated  and  experimental'’’  nonequilibrium  factors,  i.e.,  the  ratios  of  nonequilibrium 
and  equilibrium  dissociation  rates,  in  pure  nitrogen  versus  the  ratio  of  vibrational  and  translational-rotational 
temperatures. 
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•  S.O.  Macheret,  “Nonequilibrium  Processes  in  Gases  and  Plasmas  and  Their  Applications.  ”  Vanderbilt 
University,  Nashville,  TN,  April  12, 2001 
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